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CHAPTER  I:  INTRODUCTION 

In  recent  years  the  concept  of  speed  made  good  has  come  to  be 
recognized  as  one  of  the  most  meaningful  measures  of  the  off-road 
mobility  of  a  vehicle.  Speed  made  good  is  simply  the  average  speed 
that  results  when  the  straight-line  distance  between  two  points  is 
divided  by  the  time  required  to  get  from  one  point  to  the  other, 
irrespective  of  the  actual  path  taken  by  the  vehicle.  The  maximum 
speed  that  a  vehicle  can  attain  in  a  given  terrain  is  affected  by  the 
number  and  character  of  the  terrain  deterrents  it  must  overcome  or, 
if  preferable,  avoid.  Terrain  deterrents  include  soft  soils,  slopes, 
vegetation,  obstacles,  riverine  situations,  rough  terrain,  etc. 

Experience  has  shown  that  one  of  the  most  significant  single 
features  influencing  a  vehicle's  speed  is  that  of  ride  dynamics, 
that  is,  the  vibratory  activity  of  a  vehicle  caused  by  terrain  rough¬ 
ness.  Because  of  this,  a  great  deal  of  interest  has  been  shown  in  the 
ride  problem.  However,  to  date,  there  is  no  satisfactory  means  of 
relating  a  measure  of  ride  to  a  measure  of  terrain  roughness.  There 
has  actually  been  much  work  done  in  the  area  of  random  vibration 
studies,  and  the  literature  available  is  quite  abundant.  However, 
for  the  most  part,  this  work  is  largely  of  an  academic  nature  and 
oriented  toward  analysis  of  linear  systems  using  the  concepts  of 
transfer  functions.  The  results  are  very  general  and  nebulous  in 
nature,  and  an  adequate  understanding  of  them  by  an  individual  would 
require  a  substantial  background  of  education  and  experience  in 
engineering,  physics,  mathematics,  statistics,  and  electronics.  It  is 
for  this  reason  that  widespread  practical  application  of  these  results 
has  been  stymied. 
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The  General  Problem 

The  description  and  evaluation  of  vehicle  dynamic  responses  in¬ 
volves  essentially  a  study  and  interpretation  of  waveforms.  Generally, 
an  arbitrary  waveform  (the  terrain)  serves  as  an  input  to  excite  some 
deterministic  system  (a  vehicle  or  a  model) ,  and  the  outputs  are  gener¬ 
ated  as  continuous  waveforms  representing  the  resulting  motions  of  the 
system.  Obviously,  these  continuous  traces  or  waveforms  must  be  ana¬ 
lyzed  in  accordance  with  some  specified  standard  to  provide  a  basis  for 
comparisons.  This  means  simply  that  numbers  must  be  determined  that 
adequately  describe  the  pertinent  characteristics  of  each  waveform. 
Generally,  these  waveforms  are  not  well-defined  mathematical  functions, 
such  as  sinusoids,  sawtooth  waves,  etc.,  but  vary  in  seemingly  erratic 
fashions.  However,  regardless  of  how  erratic  these  waveforms  may 
appear,  their  peculiar  characteristics  can,  possibly,  be  described  ade¬ 
quately  by  statistical  methods. 

Most  engineers  are  familiar  enough  with  the  fundamentals  of  sta¬ 
tistics  to  know  that  statistical  quantities  evolve  theoretically  from 
idealized  mathematical  relations,  such  as  infinite  integrals,  that 
serve  as  approximations  to  the  real  world.  Therefore,  when  statistics 
is  applied  to  real-world  problems,  only  estimates  of  the  true  values 
can  be  obtained,  and  the  deviation  of  these  estimates  from  the  true 
values  depends  on  how  well  certain  criteria  have  been  met.  For  example, 
when  a  coin  is  flipped  a  sufficient  number  of  times,  about  as  many  heads 
as  tails  result.  This  relation  improves  as  the  number  of  flips  is 
increased,  and  theoretically  the  number  of  heads  and  tails  is  the  same 
after  an  infinite  number  of  flips. 
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The  nature  of  the  quantities  one  might  choose  to  describe  vehicle- 

terrain  responses  and  how  they  may  be  interpreted  to  define  certain 

properties  that  have  physical  meaning  in  mobility  studies  are  discussed 

here  without  rigorous  details  of  probability  and  sampling  theory. 

Averages,  measures  of 
mean,  and  dispersion 

An  average  is  a  value  that  is  typical  or  representative  of  a  set 
of  data.  Since  such  typical  values  tend  to  lie  centrally  within  a  set  of 
data  arranged  according  to  magnitude,  averages  are  also  called  measures  of 


central  tendency.  Some  of  the  most  commonly  used  averages  are  the  mean, 
the  median,  and  the  mode. 

Variation  or  dispersion  is  a  measure  of  the  degree  that  numerical 
data  tend  to  spread  about  some  reference  value.  Generally,  this  refer¬ 
ence  is  chosen  to  be  the  mean.  Several  commonly  used  measures  of  vari¬ 
ation  or  dispersion  are  the  mean  deviation,  the  standard  deviation,  the 
variance  (which  is  simply  the  square  of  the  standard  deviation) ,  the 
mean  square,  and  the  root  mean  square  (RMS).  As  an  example  of  how  such 
values  represent  data,  examine  the  four  groups  of  numbers  in  the  follow¬ 


ing  illustration: 
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All  four  groups  have  the  same  arithmetic  mean,  and  obviously,  the 
numbers  in  groups  (1)  and  (2)  are  more  uniform  than  those  corresponding 
in  (3)  and  (4).  It  is,  therefore,  apparent  that  the  arithmetic 
mean  gives  only  some  measure  of  central  tendency  and  provides  no 
indication  of  the  uniformity  of  numbers.  However,  the  numbers  in 
the  right-hand  columns  in  the  above  groups  reflect  the  squared  values 
of  the  left-hand  numbers,  and  the  mean  of  these  numbers  (the  mean 
square)  is  not  the  same  in  all  columns,  but  increases  in  magnitude 
with  an  increase  in  nonuniformity.  Hence,  the  mean  square  value 
gives  some  measure  of  nonuniformities.  Thus,  the  square  root  of 
these  mean  square  values  (RMS)  would  yield  a  similar  descriptor 
of  nonuniformity. 

Suppose  now  that  the  four  left-hand  columns  of  numbers  have  dif¬ 
ferent  mean  values.  To  provide  a  number  with  a  more  consistent  com¬ 
parison  capability,  the  deviation  about  each  mean  value  can  be  found. 
This  quantity  is  referred  to  as  the  mean  deviation.  To  illustrate  this, 
consider  again  the  four  left-hand  columns,  each  with  an  arithmetic  mean 
of  7.5.  The  mean  deviation  for  group  1  is 

I  ?•  5-7 . 5 1  +  [  7 .5-7^5  |  +  p  .  0  (1) 

6 

i.e.,  there  is  zero  dispersion  or  variation  of  the  numbers  about  the 
arithmetic  mean.  However,  for  group  4 

| 2-7. 5 1+| 13-7. 5 |+1 4-7. 5 |+|8-7. 5 1+| 10-7.51+ | 8-7.51  = 

6 

5.5  +  5.5  +  3.5  +  0.5  +  2.5  +  0.5  =  18.0 
6  6 
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(2) 


This  indicates  that  there  is  less  dispersion  in  group  1  than  in  group  4, 
as  it  should. 

A  somewhat  similar  measurement,  called  the  standard  deviation,  de¬ 
notes  the  root  mean  square  of  the  deviations  from  the  mean  and  is  some¬ 
times  called  the  RMS  deviation.  The  standard  deviation  of  the  numbers 
in  the  left-hand  column  of  group  1  is  obviously  zero  since  each  of  the 
deviations  from  the  mean  is  zero.  However,  the  standard  deviation  for 
group  4  is 


The  basic  difference  between  the  standard  deviation  and  the  mean  devia¬ 
tion  is  the  greater  influence  of  extreme  values  on  the  standard  devia¬ 
tion  than  on  the  mean  deviation.  If  the  nature  of  the  problem  is  such 
that  the  influence  of  extreme  values  are  deemed  important,  then  the 
standard  deviation  is  the  preferable  measure  of  dispersion.  The  square 
of  the  standard  deviation  is  called  the  variance.  For  normally  distrib¬ 
uted  data,  the  standard  deviation  and  variance  play  a  significant  part 
in  probability  and  statistical  estimation  theory. 

Applicability  to  terrain  roughness 

To  see  how  such  measures  might  be  used  to  describe  terrain  roughness, 
suppose  a  terrain  elevation  profile  along  a  line  is  obtained  and  plotted 
in  the  continuous  waveform  as  follows: 
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The  dispersion  measurements  used  previously  will  yield  a  number  that  de¬ 
scribes  the  deviation  of  the  elevation  amplitudes  from  the  mean  value  or 
from  the  level  datum,  if  desired.  This  quantity  should  then  provide  a 
good  measure  of  roughness.  With  this  method,  a  large  value  of  the  stan¬ 
dard  deviation  or  variance  indicates  a  large  dispersion  of  the  elevea- 
tions  and,  therefore,  a  large  amount  of  terrain  roughness;  while  a  low 
value  indicates  a  relatively  smooth  terrain.  If  the  mean  value  is  zero, 
then  the  RMS  and  standard  deviation  are  the  same,  and  the  mean  squared 
value  and  the  variance  are  the  same. 

However,  this  means  for  describing  terrain  roughness  is  not  quite 
adequate  for  evaluating  the  effect  terrain  roughness  has  on  vehicle 
response  because  a  vehicle  is  a  finely  tuned,  frequency-dependent  system 
with  a  different  response  for  each  different  frequency  of  a  sinusoidal 
displacement.  This  implies  that  a  correlation  of  a  vehicle's  response 
to  an  input  terrain  requires  a  knowledge  of  the  dispersion  of  the  ter¬ 
rain  amplitudes  per  unit  frequency.  This  will  yield  a  spectrum  of  dis- 
persons  which  when  summed  will  yield  a  single  number  representing  the 
total  dispersion.  Such  a  spectrum,  called  a  power  spectral  density 
function,  serves  the  more  flexible,  dual  purpose  of  describing  terrain 
roughness  in  a  manner  that  permits  an  evaluation  of  vehicle  response. 

A  cross-country  terrain  profile  can  be  thought  of  as  a  random  function 
of  elevation  values  versus  horizontal  distance,  i.e.  there  is  no  way  to 
predict  an  exact  elevation  at  some  given  distance.  This  randomness 
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becomes  even  more  evident  when  consideration  is  given  to  the  travel  of 
a  vehicle  across  an  off-road  area  in  which  no  specific  path  of  travel 
is  defined.  In  the  same  sense  that  a  sinusoid  can  be  defined  by  its 
amplitude  and  frequency,  a  random  profile  can  be  defined  by  its  dis¬ 
tribution  of  amplitudes  and  distribution  of  frequencies.  Representative 
measurements  have  shown  that  the  amplitudes  of  many  terrain  and  road 
profiles  can  be  reasonably  approximated  with  a  Gaussian  or  normal  dis¬ 
tribution.  ^  Although  such  a  distribution  is  not  a  stringent  require¬ 
ment  for  spectral  analysis  of  profiles,  it  enhances  the  statistical 
analysis  by  permitting  the  use  of  the  many  concepts  defined  for  normal 
distributions.  Regarding  the  distribution  of  frequencies  for  most 
profiles,  the  height  of  the  profile  is  approximately  proportional  to 
the  wavelength  of  the  frequency  being  considered.  However,  in  deter¬ 
mining  a  power  spectral  density  from  a  profile,  estimating  the  true 
function  from  a  finite  sample  and  determining  the  validity  of  the 
estimates  are  problems.  A  major  consideration  which  is  dealt  with 
later  in  the  text  is  the  Weiner-Khintchine  relation  which  relates  the 
autocorrelation  function  and  the  power  spectral  density  as  Fourier 

2 

transform  pairs,  and  is  valid  only  for  random,  stationary  profiles. 
Stationary  means  that  the  statistical  properties  of  the  profile  do  not 
change  with  position.  While  most  profiles  are  not  stationary  in  the 
strictest  sense,  valid  estimates  of  power  spectral  density  can  still  be 
made  from  profiles  that  are  quasi-stationary . 

The  general  problem  of  the  study  reported  herein  was  to  investigate 
the  validity  and  practicality  of  applying  spectral  density  methods  and 
the  inherent  statistics  to  the  problem  of  describing  terrain  roughness 
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and  vehicle  responses  in  terms  of  vehicle  speed. 

Obj  ective 

The  specific  objective  of  this  study  was  to  define  and  formulate 
the  procedures  for  implementing  a  practical,  compact  description  of 
terrain-vehicle  speed  (TVS)  systems  and  thus  establish  a  methodology 
for  relating  measures  of  ride  to  measures  of  terrain  roughness.  Using 
this  methodology,  reliable  estimates  can  be  made  of  the  average  speed 
a  given  vehicle  can  attain  over  given  terrain  as  well  as  of  its  prob¬ 
abilities  of  exceeding  certain  specified  levels  of  vibration. 

This  may  provide  a  more  practical  means  of  including  the  "ride" 
problem  in  computerized  analytical  models  used  to  quantify  the 
terrain-vehicle-driver  system,  and  thus  provide  the  military  with  a 
reliable  analytical  tool  for  (a)  developing  more  effective  vehicles 
at  less  cost  and  (b)  rating  them  on  a  competitive  basis  in  accordance 
with  the  tasks  to  be  performed. 

Scope  of  Study 

This  study  was  entirely  a  computerized  effort.  Digital  computer 
simulations  were  conducted  in  which  a  vehicle  was  run  at  several 
selected  speeds  over  terrain  profiles  with  various  levels  of  roughness. 
The  terrain  profiles  were  generated  by  computer  programs  that  con¬ 
structed  and  shaped  sequences  of  random  normal  numbers  to  provide 
representative  profiles  with  the  desired  power  spectrum,  variance,  mean 
value,  standard  deviation,  and  RMS.  These  programs  also  performed  the 
necessary  statistical  checks  for  normality,  stationarity ,  and 


randomness . 
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A  two-dimensional,  nonlinear  mathematical  model  of  an  M37  truck 
served  as  the  vehicle  in  these  simulations.  This  model  allowed  for 
the  nonlinearities  inherent  in  large  rotational  motions,  the  suspension 
elements,  bump  stops,  loss  of  ground  contact,  and  the  tire  compliance, 
which  was  represented  by  a  cluster  of  radially  projecting  springs. 

The  analysis  is  presented  chiefly  in  the  form  of  graphs  comparing 
input  and  output  statistics  and  distribution  functions.  The  input  sta¬ 
tistics  consisted  of  terrain  roughness  as  measured  by  RMS  elevation. 

The  output  is  represented  as  the  RMS  vertical  acceleration  of  the 
vehicle's  center  of  gravity. 

A  family  of  curves  was  developed  relating  vehicle  response  to 
terrain  roughness  for  each  speed  of  vehicle  traversal.  Cross  plots  of 
these  established  relations  provided  a  more  useful  relation  between 
vehicle  response  and  speed  for  various  degrees  of  terrain  roughness. 
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CHAPTER  II:  DESCRIPTION  OF  COMPUTER  CALCULATIONS 

Digital  Simulation  of  Low-Pass  Gaussian  Profiles 

An  essential  requirement  for  this  study  was  the  capability  to 
generate  artificial  terrain  profiles  composed  of  random,  normally  dis¬ 
tributed  amplitudes  with  a  given  degree  of  stationarity  and  control  of 
the  frequency  content  and  pertinent  statistics.  Specifically,  the  first 
aim  was  to  construct  a  group  of  terrain  profiles  that  were  different 
deterministically  but  were  the  same  statistically,  and  then  to  system¬ 
atically  vary  certain  statistics. 

A  series  of  computer  programs  and  subroutines,  appropriately  inter¬ 
woven,  generated  these  low-pass,  Gaussian  profiles  with  a  zero  mean  and 
specified  standard  deviation  by  the  following  procedures:  A  random- 
number  chain,  generated  by  the  power  residue  method,  was  entered  at  an 
arbitrary  starting  point,  and  12  uniform,  random  numbers  were  computed 
and  summed.  By  the  central  limit  theorem,  a  single  random,  normal 
number  was  then  computed  by  the  formula: 

V  =  (A  -  6)(an) 

where 

V  =  the  random,  normal  number 
A  =  the  sum  of  12  uniform,  random  numbers 
cr^  =  the  standard  deviation  desired  for  the  sequence  of  random, 
normal  numbers. 

This  technique  for  computing  random,  normal  numbers  is  often  referred 
to  as  the  "sum  of  the  uniform  deviates"  method. ^ 
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Generally,  the  calculations  were  repeated  until  a  sequence  of 

1000  random,  normal  numbers  was  obtained.  Although  the  mean  of  the 

resulting  sequence  was  always  very  nearly  zero,  a  shifting  operation 

was  performed  to  insure  a  mean  value  of  absolute  zero.  The  frequency 

distribution  of  the  sequence  at  this  point  approximates  that  of  white 

Gaussian  noise  in  that  its  spectral  density  is  uniformly  distributed 

2 

over  all  frequencies  with  an  average  mean  square  level  of  on  .  To 
obtain  the  desired  spectrum,  the  sequence  was  then  passed  through  a 
numerical  system  that  simulated  an  analog  low-pass  filter  with  a  cer¬ 
tain  "time  constant"  T  =  RC  and  cutoff  frequency  a  . ^  The  analog 
equivalent  of  this  system  is  shown  in  Fig.  1.  This  system  is  formu¬ 
lated  numerically  by  the  following  formula: 


i+1 


-  xi+i  (0>V  + 


Y.e 

l 


-ax 


where 

X(0,aN)  =  sequence  of  random,  normal  numbers  previously  calculated 
a  =  the  system's  cutoff  frequency  (temporal  or  spatial) 
x  =  the  time  (or  distance)  interval  between  points  in  the 
sequence 

(NOTE:  The  points  are  assumed  to  be  equidistant.) 

Y  =  the  resulting  sequence. 

The  ax  product  gives  complete  control  over  the  spectrum  shaping.  It 
was  determined  after  several  trials  that  ax  ~  0.055  gave  the  best 
normality  condition.  The  desired  density  level  is  achieved  by  speci¬ 
fying  the  appropriate  standard  deviation  a  of  the  resulting  sequence. 
This  value  is  obtained  by  the  formula: 


Figure  1.  Low  Pass  RC  Filter 
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a  =  —  — 

A  -  e~2aT 

The  computer  program  entitled  "NOISE,"  listed  in  Appendix  A,  performs 
the  operations  for  generating  the  random,  normal  sequences. 

Tests  of  Fundamental  Assumptions 

Tests  for  stationarity 

The  sequence  of  numbers  representing  a  random  terrain  profile  was 
divided  into  N  equal  intervals,  where  the  data  in  each  interval  were 
considered  to  be  independent  of  those  in  any  of  the  other  intervals. 

The  mean  and  mean  square  values  were  computed  for  each  interval.  The  re¬ 
sulting  sequences  of  the  mean  and  mean  square  values  were  tested  for 

3 

underlying  trends  by  two  tests,  the  run  test  and  the  trend  test. 
Generally,  the  run  test  is  more  powerful  for  detecting  fluctuating 
trends,  and  the  trend  test  is  more  powerful  for  detecting  underlying 
monotonic  trends.  Stationarity  required  that  the  sample  pass  both  the 
run  and  the  trend  tests  at  the  95  percent  confidence  level.  The  se¬ 
quence  of  mean  square  values  were  each  assumed  to  be  independent 

samples  of  a  random  variable  with  a  true  mean  y  and  mean  square  value 
2 

.  If  this  hypothesis  were  true,  the  variations  in  the  sequence  of 
sampled  values  would  be  random  and  display  no  trends.  Hence,  the 
number  of  runs  in  the  sequence  relative  to  some  value  (the  mean)  would 
be  as  expected  for  a  sequence  of  independent  random  observations  of  the 
random  variable,  as  given  in  Table  1.  Moreover,  the  number  of  reverse 
arrangements  in  the  sequence  would  be  as  expected  for  a  sequence  of 


14 


Table  1 

Percentage  Points  of  Run  Distribution 


a 


n  =  N/2 

0.99 

0.975 

0.95 

0.05 

0.025 

0.01 

5 

2 

2 

3 

8 

9 

9 

6 

2 

3 

3 

10 

10 

11 

7 

3 

3 

4 

11 

12 

12 

8 

4 

4 

5 

12 

13 

13 

9 

4 

5 

6 

13 

14 

15 

10 

5 

6 

6 

15 

15 

16 

11 

6 

7 

7 

16 

16 

17 

12 

7 

7 

8 

17 

18 

18 

13 

7 

8 

9 

18 

19 

20 

14 

8 

9 

10 

19 

20 

21 

15 

9 

10 

11 

20 

21 

22 

16 

10 

11 

11 

22 

22 

23 

18 

11 

12 

13 

24 

25 

26 

20 

13 

14 

15 

26 

27 

28 

25 

17 

18 

19 

32 

33 

34 

30 

21 

22 

24 

37 

39 

40 

35 

25 

27 

28 

43 

44 

46 

40 

30 

31 

33 

48 

50 

51 

45 

34 

36 

37 

54 

55 

57 

50 

38 

40 

42 

59 

61 

63 

55 

43 

45 

46 

65 

66 

68 

60 

47 

49 

51 

70 

72 

74 

65 

52 

54 

56 

75 

77 

79 

70 

56 

58 

60 

81 

83 

85 

75 

61 

63 

65 

86 

88 

90 

80 

65 

68 

70 

91 

93 

96 

85 

70 

72 

74 

97 

99 

101 

90 

74 

77 

79 

102 

104 

107 

95 

79 

82 

84 

107 

109 

112 

100 

84 

86 

88 

113 

115 

117 
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independent  random  observations  of  the  same  variable,  as  given  in 
Table  2.  If  the  number  of  runs  or  reverse  arrangements  was  signifi¬ 
cantly  different  from  the  expected  number  given  in  the  tables,  the 
hypothesis  of  stationarity  would  be  rejected.  The  details  of  these  two 
tests  are  presented  here  to  give  an  understanding  of  the  principles 
involved  in  the  calculations. 

Run  test.  Consider  the  sequences  of  mean  values  X_^  and  mean 

square  values  Y^  for  a  given  profile  to  be  classified  into  mutually 

exclusive  categories  identified  by  plus  (+)  or  minus  (-) ,  where 

X.  >X=+  ,  X.<  X  =  -  ;  and  likewise  for  the  Y.  .  A  run  is  de- 
11  l 

fined  as  a  sequence  of  identical  observations.  The  number  of  runs 
that  occur  in  a  sequence  of  observations  indicates  whether  or  not 
results  are  independent  random  observations. 

Let  it  be  hypothesized  that  there  is  no  trend  in  the  sequences, 
that  is,  no  variations  in  the  mean  and  mean  square  values,  which  implies 
stationarity.  If  the  sequences  of  mean  and  mean  square  values  are  inde¬ 
pendent  observations  and  the  number  of  (+)  observations  equals  the  number 
of  (-)  observations,  then  the  number  of  runs  in  the  sequence  will  have  a 
sampling  distribution  as  given  in  Table  1.  The  hypothesis  can  be  tested 
to  any  desired  level  of  significance  a  by  comparing  the  observed  runs 

in  the  interval  between  r  n  and  r  ,  where  n  =  N/2  ,  N 

n;  1  -  a/2  n;  a/2 

being  the  total  number  of  observations.  If  the  observed  runs  fall  out¬ 
side  the  interval,  the  hypothesis  is  rejected  at  the  a  level  of 
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Table  2 

Percentage  Points  of  Reverse  Arrangement  Distribution 


a 


N 

0.99 

0.975 

0.95 

0.05 

0.025 

0.01 

10 

9 

11 

13 

31 

33 

35 

12 

16 

18 

21 

44 

47 

49 

14 

24 

27 

30 

60 

63 

66 

16 

34 

38 

41 

78 

81 

85 

18 

45 

50 

54 

98 

102 

107 

20 

59 

64 

69 

120 

125 

130 

30 

152 

162 

171 

263 

272 

282 

40 

290 

305 

319 

460 

474 

435 

50 

473 

495 

514 

710 

729 

751 

60 

702 

731 

756 

1013 

1038 

1064 

70 

977 

1014 

1045 

1369 

1400 

1437 

80 

1299 

1344 

1382 

1777 

1815 

1360 

90 

1668 

1721 

1766 

2238 

2283 

2336 

100 

2083 

2145 

2198 

2751 

2804 

2865 
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significance.  For  an  example,  consider  the  sequence  of  the  observa¬ 
tions  below: 


(1) 

5.2 

(6) 

5.7 

(ID 

6.2 

(16) 

5.4 

(2) 

5.1 

(7) 

5.0 

(12) 

6.4 

(17) 

6.4 

(3) 

5.7 

(8) 

6.5 

(13) 

4.7 

(18) 

5.8 

(4) 

5.2 

(9) 

5.4 

(14) 

5.4 

(19) 

6.7 

(5) 

4.3 

(10) 

5.8 

(15) 

5.9 

(20) 

5.2 

The  mean  value  of  the  observations  is  5.6.  Let  all  observations  with  a 


value  greater  than  5.6  be  (+)  and  those  less  than  5.6  be  (-) .  The 
result  is 


4* 


i 


2 


5 


+ 

6 


~ 

+++ 

— 

+ 

+++ 

- 

7 

8 

9 

10 

11  12 

13 

There  are  13  runs  and  20  observations.  Let  it  be  hypothesized  that  the 
observations  are  independent.  The  acceptance  region  for  this  hypothesis 
is 


r 10 ;  1  -  a/2 


<  r  <  r 


10;  a/2 


From  Table  1  for  a  =  0.05,  r^.  x  _  a/2  =  r^.  9?5  =  6  ,  and  r^. 

=  riri  .  noc  =  15  .  The  hypothesis  is  accepted  because  13(r)  falls 
0t/^-  XvJj  \J  •  U  ^  -3 

between  6  and  15.  Therefore,  the  run  test  is  passed  and  there  is  no 
evidence  of  an  underlying  fluctuating  trend. 

Trend  test.  The  trend  test  counts  the  number  of  times  that  x.  > 


-  l 

Xj  for  i  <  j  .  Each  such  inequality  is  called  a  reverse  arrangement. 
Let  the  total  number  of  reverse  arrangements  be  denoted  by  A  .  Test 
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the  same  sequence  of  20  observations  used  in  the  previous  example  for  a 
trend  at  the  a  =  0.05  level  of  significance.  For  example,  the  first 
number,  5.2,  is  greater  than  four  of  the  subsequent  numbers.  The  number 
of  reverse  arrangements  would  therefore  be  4.  The  second  number,  5.1, 
is  greater  than  three  of  the  subsequent  numbers,  with  the  first  number 
excluded  from  this  comparison.  The  number  of  reverse  arrangements  is  3. 
The  third  number,  5.7,  is  greater  than  eight  of  the  subsequent  numbers, 
with  the  first  two  numbers  excluded.  Thus,  the  number  of  reverse 
arrangements  is  as  follows: 


A1 

=  4 

A  =  6 

6 

> 

H* 

II 

6 

A16 

A2 

=  3 

A?  =  1 

A12 

6 

A17 

A3 

=  8 

A8  "  10 

A13  = 

0 

A18 

\ 

=  3 

A9  =  2 

> 

I—1 

II 

1 

A19 

A5 

=  0 

aio  =  4 

A15  = 

3 

The  total  number  of  reverse  arrangements  is  62.  Now  let  it  be 
hypothesized  that  the  observations  are  independent  observations  where 
there  is  no  trend.  The  acceptance  region  for  this  hypothesis  is 

A20;  1  -  a/2  *  A  <  A20;  a/2 


From  Table  2  for  a  =  0.05,  A  n  ,  =  A  =  64  and 

20;  1  -  a/2  20;  0.975 

A20‘  a/2  =  A20‘  0  025  =  ^5  *  Hence>  since  A  =  62  falls  outside  this 
interval,  the  hypothesis  is  rejected  at  the  5  percent  level  of  signifi¬ 
cance.  A  plot  of  the  20  observations  (Fig.  2)  shows  a  slight  monotonic 


increasing  trend  that  was  not  detected  by  the  run  test.  The  straight 
line  delineating  the  trend  was  determined  by  a  least  squares  fit  of  the 
data.  The  sensitivity  of  the  test  at  this  level  of  significance  is 


20 

perceived  from  the  small  value  of  the  slope  of  the  regression  line. 
Therefore,  for  a  sequence  such  as  this,  the  hypothesis  of  stationarity 
would  be  rejected. 

Test  for  randomness 

Randomness,  the  most  fundamental  concept  underlying  all  the  theory 
of  estimation  and  testing  of  hypotheses,  is  far  from  easy  to  understand. 
It  is  generally  interpreted  as  a  lack  of  bias,  implying  no  preference 
for  any  single  element  or  group  of  elements.  The  determination  of  this 
lack  of  bias  is  centered  around  certain  specific  elements  of  bias  (non¬ 
randomness)  that  can  be  detected  from  an  analysis  of  the  internal 
structure  of  the  sample.  This  can  be  accomplished  by  recording  the 
individual  measurements  that  make  up  a  sample  in  the  precise  order  in 
which  they  were  originally  obtained  and  performing  the  run  test  used 
previously  to  test  stationarity.  However,  this  method  is  not  practical 
(or  necessary)  for  samples  containing  a  large  number  of  points.  Having 
passed  the  test  for  stationarity,  a  test  for  randomness  effectively 

reduces  to  a  test  for  the  presence  of  sinusoids  due  to  periodic  or 

4 

almost  periodic  components  in  the  data.  The  most  powerful  method  of 
detecting  sinusoidal  components  in  otherwise  random  data  is  presented 
by  an  autocorrelogram,  i.e.  a  plot  of  the  autocorrelation  function. 

For  any  purely  random  data,  the  autocorrelation  function  will  always 
approach  a  value  equal  to  the  square  of  the  mean  value  as  the  lag 
becomes  large.  On  the  other  hand,  the  autocorrelation  function  for  a 
sine  wave,  or  collection  of  sine  waves,  will  continue  to  oscillate  no 
matter  how  large  the  lag  becomes.  Details  concerning  the  autocorre¬ 
lation  function  and  its  computation  are  presented  later  in  this  report. 
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Test  for  normality 

Normality  is  tested  conveniently  by  the  "Chi-square  Goodness-of- 
the-Fit-Test . "  The  general  procedure  of  the  test  involves  the  use  of  a 
statistic  with  an  approximate  Chi-square  distribution  as  a  measure  of 
the  discrepancy  between  an  observed  probability  density  function  and  the 
theoretical  probability  density  function.  The  approximate  Chi-square 
value  of  the  given  data  is  obtained  by: 


2 

X 


K 


1 

i=l 


V 


2 


F. 

l 


Where  f  is  the  frequency  of  occurrence  in  the  subgroup  interval  of 

the  original  sequence,  K  is  the  total  number  of  subgroups  into  which 

the  sequence  is  divided,  and  F^  is  the  number  of  expected  occurrences 

of  the  theoretically  normal  distribution.  The  Chi-square  goodness-of- 

fit  is  influenced  by  the  choice  of  the  number  of  subgroups,  K  .  For  the 

case  where  the  test  will  be  performed  at  a  =  0.05  level  of  signifi- 
9 

cance,  William  suggests  that  the  minimum  number  of  subgroups  K  be 
determined  by  the  formula: 

K  =  1 . 8 7 (N  -  1)2/5  (4) 

where  N  is  the  total  number  of  points  in  the  sequence  (1000  for  this 
study) .  To  visualize  the  magnitude  of  K  for  a  given  N  ,  a  limited 


tabulation  is  as  follows: 
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N 


100 

200 

300 

400 

500 

600 

700 

800 

900 

1000 


1500 


2000 

The  region  of  acceptance  at  the  a 


K_ 

12 

16 

18 

21 

22 

24 

26 

27 

28 
30 

I 

I 

I 

35 

I 

I 

I 

39 

0.05  level  of  significance  is 


2 

X 


< 


2 


X 


n:a 


2 

where  X  n.a  is  the  Chi-square  value  of  the  theoretically  normal  dis¬ 
tribution  function  at  the  significance  level  a  . 

A  computer  subroutine  program  was  developed  to  test  the  normality 
of  sequences  of  any  length  by  means  of  this  Chi-square  concept.  This 
program  is  a  part  of  the  STANOR  program  (Appendix  A)  that  performs  all 
the  tests  of  fundamental  assumptions  previously  mentioned. 


Autocorrelation  Function 


Generally,  the  principal  application  for  an  autocorrelation 
function  measurement  of  physical  data  is  to  establish  the  dependence  of 
values  upon  one  another.  However,  the  main  interest  here  is  that  the 
autocorrelation  function  affords  a  means  of  computing  the  power  spectral 
density  function.  The  autocorrelation  function  is  found  by  taking  two 
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identical  waveforms,  repeatedly  shifting  one  with  respect  to  the  other, 
and  then  calculating  the  correlation  coefficient  for  each  shift.  The 
correlation  coefficient  for  a  given  shift  or  lag  is  obtained  by  multi¬ 
plying  the  ordinates  of  the  two  waveforms  at  corresponding  points  on  the 
horizontal  scale  and  averaging  the  products  over  a  sufficiently  long 
range,  usually  several  periods  of  the  lowest  frequency  present,  to  pro¬ 
vide  a  steady  average.  This  is  illustrated  in  Fig.  3a.  The  products 
a^b^  ,  & 2^ 2  ’  a3^3  *  etc‘»  are  computed,  added  together,  then  divided 
by  the  total  number  of  products  in  the  sample.  The  correlation  function 
is  obtained  from  a  plot  of  the  correlation  coefficient  as  a  function  of 
the  shift  or  lagged  values  (Fig.  3b).  The  shift  values  may  be  in  terms 
of  time  or  distance,  depending  upon  whether  the  original  waveforms  are 
functions  of  time  or  distance.  The  autocorrelation  function  has  a  maxi¬ 
mum  value  for  zero  shift,  which  is  equal  to  the  mean  square  value  of  the 
waveform,  and  decreases  to  zero,  if  the  mean  value  of  the  waveform  is 
assumed  to  be  zero,  for  infinite  value  of  the  shift  unless  the  waveform 
contains  periodic  components.  If  the  waveform  has  periodic  components, 
the  autocorrelation  functions  also  has  periodic  components  of  the  same 
period.  The  autocorrelation  functions  for  several  common  waveforms  are 
shown  in  Fig.  3c. 

The  autocorrelation  function  of  a  discrete  sequence  is  defined  by 

Vx)  =  n™  n  x(k)  *  x(k  +  t)  (5) 

where  R  (x)  is  the  value  of  the  autocorrelation  function  when  the 

X 

lagged  distance  is  x  and  xOO  is  the  original  sequence.  The  dis¬ 
crete  form  of  equation  5  is  a  direct  translation  of  the  autocorrelation 


Figure  3a.  Calculation  of  Correlation  Coefficient 

At  corresponding  points  on  horizontal  scale,  the  ordinates  of  the 
two  waveforms  are  multiplied.  The  correlation  coefficient  is 
the  average  of  the  a^b^  products. 


Figure  3b,  The  Correlation  Function 
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function  for  a  continuous  and  infinite  sequence  given  by  equation  6: 


VT> 


lim  1 
T-**>  x 


fT 
■  0 


X (t)  X (t  +  x)dt 


(6) 


This  equation  is  true  for  a  sequence  that  has  an  infinite  length.  How¬ 
ever,  for  sequences  with  finite  length,  consistent  results  require  that 
the  autocorrelation  function  (equation  5)  be  modified  to  the  form 

1  N-t 

r(t)  =  rr—  l  xOO  ‘  X(k  +  x)  (7) 

This  equation  differs  from  equation  5  in  that  the  quantity  x(k)  * 

X'(k  +  t)  is  averaged  by  dividing  by  the  actual  length  over  which  x(k) 
and  x(k  +  x)  are  multiplied. 

For  a  random  sample  with  a  zero  mean,  this  function  should  approach 
zero  and  remain  at  zero  as  the  lag  values  increase,  and  the  function 
should  not  be  periodic.  (See  Fig.  3c  for  random  wave.)  The  extent  to 
which  the  function  fulfills  these  conditions  indicates  whether  the  wave¬ 
form  (terrain  profile)  is  sufficiently  random  to  be  described  by  its 
power  spectral  density.  However,  an  unsatisfactory  autocorrelation 
function  can  frequently  be  improved  by  performing  a  filtering  operation 
(detrending)  on  the  original  data  to  remove  the  influence  of  longer 
wavelengths.  A  most  useful  and  very  significant  property  of  the  auto¬ 
correlation  function  is  that  it  happens  to  be  the  Fourier  transform  of 
the  power  density  spectrum  of  a  random,  stationary  waveform.  Therefore, 
having  once  obtained  the  autocorrelation  function,  it  is  then  possible 
to  compute  the  power  spectral  density  function. 
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Power  Spectral  Density 

The  power  spectral  density  (PSD)  is  a  second  moment  or  variance 
density  spectrum.  PSD  of  a  random,  stationary  function  is  defined  as 
the  Fourier  transform  of  its  autocorrelation  function  as  follows: 


P(fi)  =  2 
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where  R  (x)  is  the  value  of  the  correlation  function  for  shift  x 
x 

and  is  the  frequency  in  eye/ in.  if  x  is  measured  in  inches  and 

cyc/sec  if  measured  in  seconds.  The  integral  over  the  frequency  range 
of  the  PSD  is  the  variance.*  The  raw  estimate  of  the  power  spec¬ 
trum  is  obtained  by  the  following  numerical  method: 
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=  autocorrelation  value  of  series  x  (the  waveform)  at 
lag  x 

=  raw  power  spectral  estimate  of  series  x  at  frequency 
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*  This  is  true  if  the  function  has  a  zero  mean.  It  was  shown  pre¬ 
viously  that  for  functions  with  a  zero  mean,  the  mean  square  and 
variance  are  the  same. 
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NOTE:  P^  =  P 
x  x 

At  =  constant  sampling  interval  (may  be  time  or  distance) 
m  =  maximum  lag  used  in  autocorrelation  function  not  to  exceed 
199. 

The  raw  estimates  are  then  smoothed  by  using  "Hamming"  coefficients  as 
follows : 
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where  SP^  is  the  smoothed  power  spectral  estimate  of  series  x 
llTT 

at  frequency  .  Smoothing  is  necessary  because  the  raw  estimate 

is  an  inefficient  estimate  of  the  true  spectral  density  because  the 
finite  record  length  precludes  identifying  frequencies  exactly.  A 
more  reliable  statistical  estimate  can  be  obtained  by  averaging  over 
neighboring  frequencies.  Trial-and-error  methods  have  proved  more 
powerful  in  obtaining  sets  of  averaging  coefficients  that  produce  re¬ 
liable  spectral  estimates  than  has  any  particular  theory.  There  are 
several  well-established  sets  of  coefficients  available,  and  in  general, 
the  Hamming  method  is  quite  satisfactory.  The  computer  program  "PSD" 
(Appendix  A)  computes  both  the  autocorrelation  and  the  PSD  estimates. 

A  plot  of  two  terrain  power  spectra  is  shown  in  Fig.  4.  The 

2 

ordinates  are  in  terms  of  ft  / cyc/ft,  and  the  abscissa  are  in  units  of 
cyc/ft,  the  reciprocal  of  wavelength.  These  spectra  measure  the  rate 
of  change  of  the  mean  squared  elevation  at  each  frequency  with  fre¬ 
quency  and  are  primarily  geometric  quantities  that  do  not  contain  any 


Figure  4.  Examples  of  Terrain  Power  Spectra 
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units  of  time.  The  area  under  each  power  spectrum  represents  the  mean 
square  value  of  the  total  roughness  for  the  respective  profile.  The 
area  under  the  curve  between  any  two  selected  wavelengths  will  indicate 
the  contribution  to  the  total  roughness  that  this  range  of  wavelengths 
produces.  Thus,  from  the  terrain  power  spectrum  curve,  the  wavelengths 
that  make  significant  contributions  to  the  variation  in  the  elevation 
measurements  can  be  determined. 

The  terrain  elevation  power  spectrum  can  be  used  also  with  certain 

vehicle  characteristics  to  determine  wavelengths  and  frequencies  that 

induce  significant  forces  and  motions  if  the  velocity  of  the  vehicle  is 

introduced.  This  is  done  by  multiplying  the  abscissa  of  the  power 

spectrum  curve  by  the  vehicle  velocity  to  obtain  the  frequency  cycles 

per  second  (cps) ,  and  by  dividing  the  ordinates  of  the  power  spectrum 

2 

by  the  vehicle  velocity  to  obtain  the  units  of  ft  /cyc/sec.  A  different 
curve  obviously  results  for  each  velocity  because  changes  in' velocity 
tend  to  stretch  or  shorten  the  curve.  Now,  take  a  hypothetical  case 
where  a  vehicle  is  driven  onto  a  sinusoidal  exciter  platform  with  a 
controlled  frequency,  adjustable  amplitude,  and  electronic  transducers 
to  measure  the  resulting  reaction  forces  and  motions.  By  varying  only 
the  frequency,  a  response  of  the  force/unit  displacement  versus  fre¬ 
quency  (Fig.  5)  can  be  obtained.  Note  the  two  predominant  frequencies 
at  about  2  and  15  cps.  These  are  the  basic  natural  frequencies  of  the 
center  of  gravity  and  the  axles.  Variations  in  the  terrain  profile 
that  will  excite  frequencies  of  2  and  15  cps  will  thus  generate  large 
forces  and  motions  between  the  vehicle  and  the  terrain.  These  fre¬ 
quencies  are,  therefore,  of  interest  in  determining  the  influence  of 
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the  terrains  on  the  reaction  of  the  vehicle.  This  relation  between 

forces  (or  motion  if  so  desired)  and  frequency  can  now  be  combined  with 

the  terrain  power  spectrum  to  obtain  an  estimate  of  the  dynamic  force  or 

motion  that  the  vehicle  will  exert  on  the  terrain,  and  vice  versa.  To 

do  this,  the  terrain  spectrum  must  first  be  modified  by  introducing  the 

vehicle  velocity,  as  previously  described.  If  this  is  done,  then  to 

determine  a  power  spectrum  of  the  dynamic  force  that  the  vehicle  will 

exert  on  the  terrain  is  simply  a  matter  of  multiplying  the  two  spectra. 

It  is  worth  mentioning  that  the  vehicle  characteristics  (Fig.  5)  must  be 

squared  when  this  mathematical  operation  is  performed  to  give  the  re- 

2 

suiting  units  of  lb  /cyc/sec  in  the  ordinates  of  the  response  spectrum. 

A  typical  result  of  these  computations  can  be  perceived  from  the  plots 
in  Figs.  6  and  7  for  the  44  ft/sec  speed.  The  influence  of  the 
vehicle's  natural  frequencies  (2  and  15  cps)  is  vividly  portrayed  by 
the  dynamic  response  spectrum  (Fig.  7).  This  response  spectrum  also 
has  the  same  characteristics  as  previous  spectra,  i.e.  the  area  under 
the  curve  represents  the  mean  square  value  of  the  response.  With  this 
value,  it  is  easy  to  determine  the  root  mean  square  (effective)  value. 
This  RMS  response  value  is  the  effective  force  (or  motion)  that  the 
terrain  exerts  on  the  vehicle  or  vice  versa,  for  that  velocity.  Such  a 
spectral  analysis  provides  a  great  deal  of  information  on  and  insight 
into  just  how  the  frequencies  of  the  terrain  influence  a  vehicle's 
response  and  further  pinpoints  critical  speeds  and  terrain  wavelengths. 
However,  this  analysis  is  based  on  the  fact  that  the  vehicle  is  a  linear 
system;  that  is,  if  the  amplitude  of  the  sinusoidal  shaker  is  doubled, 
the  only  effect  on  the  vehicle's  response  is  the  doubling  of  the  output 
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for  all  frequencies.  This,  however,  is  not  the  case  for  actual  vehi¬ 
cles  which,  unfortunately,  are  nonliner. 

Vehicle  Model 

A  two-dimensional,  4-deg-of-freedom  mathematical  model  for  the  M37, 
3/4-ton  truck  was  developed  for  this  study.  The  truck  was  modeled  in 
terms  of  a  series  of  coupled  second-order  ordinary  differential  equa¬ 
tions  that  are  solved  simultaneously  to  describe  the  motions  of  the 
truck  as  it  traverses  the  selected  terrain  profiles.  The  frame  of  the 
truck  was  considered  rigid,  and  only  the  pneumatic  tires  and  suspensions 
were  considered  to  contribute  to  the  sprung  motion  of  the  frame.  The 
model  included  all  the  pertinent  nonlinearities  in  the  suspensions, 
those  inherent  in  large  rotational  motions  and  loss  of  ground  contact. 
Also,  the  segmented  wheel  concept^  was  employed  to  describe  tire-terrain 
compliance.  This  added  still  another  nonlinear  feature  to  the  model. 

The  equations  and  a  schematic  diagram  for  the  model  are  as  follows: 


Schematic  of  Truck  Model 
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a.  Forces  on  body. 

Mz  =  K^(A^)(z^  -  z  -  a  sin  0)  +  C^(A^)(z^  -  z  -  a  0  cos  0) 

+  K2(A2)(z2  -  z  +  b  sin  0)  +  C^k^)(.z2  -  z  +  b  6 

cos  0)  -  Mg 

10  =  K^(A^)a(z^  -  z  -  a  sin  0)  +  C^(A^)a(z^  -  z  -  a  0  cos  0) 

-  K2(A2)b(z2  -  z  +  b  sin  6)  -  C2(A2)b(z2  -  z  +  b  0 

cos  0) 

b .  Forces  on  front  axle. 

M^z^  =  -  K^(A^) (z^  -  z  -  a  sin  0)  -  C^CA^Hz^  -  z  -  a  0 

10 

cos  0)  +  l  Y11(p11  -  z1)  -  M1g 
c ■  Forces  on  rear  axle. 

M2^2  =  ~  K2^A2^z2  “  z  +  b  sin  “  C2(A2)(z2  -  z  +  b  0 

10 

cos  9)  +  l  Yi2(Pi2  "  z2)  ~  M2g 

where 

A^  =  z^  -  z  -  a  sin  0  A^=z^-z-a0  cos  0 

A2  =  z2  -  z  +  b  sin  0  A2  =  z2  -  z  +  b  0  cos  0 

For  this  study,  the  suspension  spring  coefficients  were  presented  by 
third-order  polynomials,  as  shown  in  Fig.  8.  These  polynomials  were 
obtained  by  curve  fitting  the  actual  force-deflection  relations,  and 
they  are  reasonable  approximations.  The  suspension  damping  is  as 


follows : 
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C.,  (A  )  =  11.8  (compression) 

11  m . 

=  22.8  (extension) 

in. 

=  12.0  — — — sec  (compress-[on) 

=  49.0  t - (extension) 

in. 

No  damping  was  incorporated  in  the  tire  compliance  since  in  actual  ve¬ 
hicles  this  damping  is  negligible  compared  to  that  of  the  suspensions. 
This  truck  model  was  forced  to  traverse  each  terrain  condition  at  a 
constant  speed,  and  the  outputs  consisted  of  motions  of  the  main  frame 
and  axles  in  terms  of  displacement-,  velocity-,  and  acceleration-time 
histories.  The  differential  equations  were  programmed  for  solution  on 
the  GE-400  computer  by  the  Runge-Kutta-Gill  numerical  integration 
scheme  (Appendix  A) . 


CHAPTER  III:  TEST  PROCEDURES  AND  DATA  ANALYSIS 


39 


Similarity  of  Input  and  Output  Statistics 

The  response  of  linear  systems  to  random  excitations  is  well 
defined  in  the  form  of  mathematical  solutions  employing  the  concept  of 
transfer  functions.  For  example,  if  a  random  input  forcing  function 
can  be  assumed  to  have  a  normal  probability  distribution  about  a  mean 
value  of  zero,  then  the  response  of  this  input  will  also  be  a  random 
quantity  with  a  normal  probability  distribution  about  a  mean  value 
of  zero.  In  other  words,  the  input  statistics  (the  terrain)  are  trans¬ 
formed  in  a  linear  fashion  by  the  transfer  characteristics  of  the 
system  (the  vehicle)  to  produce  response  statistics  that  are  related  by 
the  expression 

-  M“>l2 

where 

P (to)  ,P(w) ^  =  the  output  and  input  spectral  density  functions, 

respectively 

H(w)  =  the  system  transfer  function. 

Hence,  for  a  given  input  spectral  density  function,  the  output  spectral 
density  function  is  defined  as  the  product  of  this  spectral  density 
function  and  the  square  of  the  system's  transfer  function. 

Unfortunately,  the  topic  of  nonlinear  response  to  random  excita¬ 
tions  does  not  readily  lend  itself  to  analytical  treatment.  Consequen¬ 
tly,  the  most  productive  investigations  in  this  area  have  been  empirical 
in  nature.  Very  little  is  known  about  how  random  inputs  are  transformed 
by  nonlinear,  multi-degree-of-freedom  systems. 
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Before  any  effort  was  made  to  develop  a  method  for  defining 
terrain-vehicle-speed  relations  in  a  statistical  fashion,  an  important 
question  had  to  be  answered:  "Will  constant  statistical  inputs  to  a 
nonlinear  system  produce  constant  statistical  outputs?"  A  negative 
answer  to  this  question  is  equivalent  to  saying  that  there  is  little 
hope  for  assembling  a  rational  scheme  for  statistical  analysis  of  non¬ 
linear  systems.  To  answer  this  question,  three  terrain  profiles  were 
generated  (Fig.  9)  that  were  different  deterministically,  but  were  the 
same  statistically.  That  is,  they  had  essentially  the  same  PSD f  dis¬ 
tribution  function  and  identical  RMS  values,  although  the  elevation 
sequences  composing  their  waveforms  were  quite  different.  The  con¬ 
struction  of  these  profiles  was  influenced  by  the  fact  that  a  number  of 
investigators  have  observed  that  the  power  spectra  (PSD)  of  natural 
surfaces  and  most  man-made  surfaces,  such  as  runways  and  highways,  can 
be  expressed  by  the  equation: 

P(ft)  =  Cfi“n  (9) 

Published  data  further  indicate  that  the  exponent  n  for  both 

g 

natural  and  man-made  surfaces  is  approximately  two.  Therefore,  a 
requirement  for  constructing  these  profiles  was  that  their  power  spectra 
would  be  approximated  by  equation  9  with  an  exponent  of  two.  Although 
the  nearly  20-to-l  ratio  between  the  vertical  and  horizontal  scales 
produces  distorted  plots,  the  differences  in  features  of  the  profiles 
are  apparent.  The  RMS  elevation  was  calculated  for  each  profile,  and 
in  each  case,  a  value  of  8.5  was  obtained.  The  RMS-distance  plots  are 
shown  in  Fig.  10.  These  profiles  were  purposely  made  extremely  rough 
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to  assure  that  all  the  vehicle's  nonlinearities  would  be  influenced  by 

the  maximum  degree.  Logarithmic  plots  of  the  profile  spectra  are  shown 

_2 

in  Fig.  11.  The  straight  line,  whose  equation  is  P(fi)  =  0.0785^  , 

drawn  through  each  plot  is  to  enhance  a  visual  comparison  of  the  simi¬ 
larity  among  the  three  spectra  and  illustrate  the  capability  available 
for  approximating  given  spectral  relations. 

The  truck  model  traversed  each  of  the  three  profiles  at  a  speed  of 
5  mph  (88  in. /sec).  The  vertical  acceleration-time  histories  of  the 
vehicle's  CG  and  front  axle  are  shown  in  Figs.  12  and  13,  respectively. 
Their  respective  power  spectra  are  given  in  Figs.  14  and  15.  The  corre¬ 
sponding  acceleration-time  histories  appear  similar  in  form,  but  it  is 
obvious  from  the  distribution  of  the  peaks  that  the  waveforms  are  not 
truly  stationary.  However,  if  a  quasi-stationary  condition  is  consid¬ 
ered  for  each  case,  the  response  spectra  can  be  computed.  Close  exami¬ 
nation  of  the  corresponding  response  spectra  reveal  similar  frequency 
distributions,  but  slight  deviations  in  the  low  frequency  range  are 
noted  for  the  CG  motion.  Unfortunately,  at  this  speed  this  low  fre¬ 
quency  range  corresponds  closely  to  the  basic  natural  frequency  of  the 
vehicle  CG.  Also,  the  axle  response  spectra  for  the  test  on  terrain 
No.  3  is  noticeably  lower  than  that  on  the  other  two  spectra.  These 
deviations  in  the  output  spectra  are  reflected  in  the  RMS  acceleration- 
distance  histories  (Figs.  16  and  17).  In  Fig.  16,  which  shows  the  RMS- 
distance  history  of  the  CG  motion,  the  RMS  traces  tend  to  converge  and 
reach  values  of  1.02,  0.95,  and  0.91  g,  respectively,  after  about 
4000  in.  of  travel.  The  RMS  deviations  at  this  point  are  really  quite 
small,  and  indications  are  that  they  would  converge  to  a  common  value 


;ure  11.  Logarithmic  Plots  of  the  Profile  Spectra  for  the  Computer  Generated  Profiles  with 

Similar  Statistical  Content 


Figure  12.  Acceleration^-Time  Histories  of  Vehicle's 
Center  of  Gravity  Traversing  Three  Profiles 
with  Same  Statistical  Content  at  5  mph 


Figure  13.  Acceleration-Time  Histories  of  Vehicle's 
Front  Axle  Traversing  Three  Profiles  with 
Same  Statistical  Content  at  5  mph 


Figure  16.  RMS  Acceleration  versus  Distance  of  CG  Response 
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if  the  test  were  evaluated  over  a  longer  distance. 

A  somewhat  similar  occurrence  is  noted  for  the  axle  motions 
(Fig.  17),  where  terminal  RMS  values  of  3.55,  3.38,  and  3.06  g  were 
recorded.  As  indicated  by  the  response  spectra,  the  RMS  acceleration 

at  the  axle  for  the  test  on  terrain  No.  3  was  lower  than  that  recorded 

( 

for  the  other  two.  The  noticeable  deviations  in  the  low  frequency  range 
of  the  CG  response  are  believed  to  further  indicate  an  insufficient 
length  of  sampling.  Each  of  these  tests  represented  only  about  330  ft 
of  travel  for  a  vehicle  with  a  wheel  base  of  over  9  ft. 

The  autocorrelation  functions  of  the  CG  motions  were  plotted  for 
each  of  the  three  tests,  and  the  results  are  shown  in  Fig.  18.  In 
each  test,  a  noticeable  periodicity  occurs.  This  is  believed  due 
largely  to  the  periodic  influence  of  the  front  and  rear  axles  hitting 
the  same  profile  elevations  with  a  time  delay  dependent  on  the  wheel 
base  and  the  vehicle  speed.  This  is  probably  the  chief  contributor  to 
the  nonstationarity  of  the  output.  Therefore,  in  all  probability  if 
the  test  had  been  conducted  long  enough  to  allow  these  periodic  in¬ 
fluences  to  achieve  stationarity,  it  is  conceivable  that  at  this  point 
and  beyond  the  output  RMS  results  would  coincide.  If  this  were  the 
case  and  since  this  periodic  phenomenon  is  a  function  of  the  wheel 
base,  it  appears  logical  that  a  minimum  distance  criterion  could  be 
established  in  terms  of  the  vehicle's  wheel  base  for  evaluating  vehicle 
responses  in  cross-country  environments.  For  example,  one  might  say 
that  to  statistically  evaluate  the  dynamic  response  of  cross-country 
vehicles  requires  travel  over  a  profile  that  is  suitably  stationary 
and  whose  length  is  at  least  n  times  that  of  the  wheel  base  of  the 
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Figure  17.  RMS  Acceleration  Versus  Distance  for  Front  Axle  Motion 
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Figure  18.  Autocorrelation  Functions  of  Vehicle's  CG  Acceleration 


52 


vehicle.  The  value  for  n  would  have  to  be  established  from  either 
experimental  data  or  simulations. 

Development  of  TVS  Relations 

It  was  concluded  from  the  previous  tests  that  the  vehicle  response 
statistics  obtained  from  traversing  terrains  with  constant  statistics 
and  a  given  terrain  power  spectra  will  be  similar.  This  fact  allows 
progression  to  the  next  step,  which  is  the  main  objective  in  this  study: 
to  determine  the  effect  of  terrain  roughness  and  vehicle  speed  on  the 
vehicle's  statistical  response.  Simulations  with  the  truck  model  were 
conducted  over  eight  profiles  of  varying  levels  of  terrain  roughness  at 
speeds  of  5,  7.5,  10,  15,  and  20  mph  over  each  of  the  eight  profiles. 
Only  the  RMS  elevation  of  the  profile  was  varied,  while  the  shape  of 
the  spectra  was  maintained  constant.  The  nature  of  this  variation  is 
illustrated  in  Fig.  19,  which  shows  logarithmic  plots  of  the  terrain 
profile  spectra  for  three  levels  of  roughness.  The  RMS  elevations  for 
these  three  profiles  were  2.8,  8.5,  and  63.7,  respectively.  The 
profile  with  RMS  =  63.7  would  preclude  any  type  of  vehicular  traffic 
and  is  shown  only  to  illustrate  that,  the  profile  roughness  can  be 
changed  without  altering  the  profile's  frequency  distribution.  This 
method  of  changing  only  the  profile  RMS  simply  alters  each  elevation 
point  by  a  proportionate  amount.  That  is,  if  the  original  profile  had 
an  RMS  of  5  and  the  RMS  was  changed  to  10,  the  features  of  the  resulting 
profile  would  be  identical  to  those  of  the  original,  but  the  elevations 
would  be  twice  those  of  the  original. 

The  straight  lines  drawn  through  the  three  spectra  each  represent 
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Figure  19.  Power  Spectra  for  Computer  Generated  Terrains 
of  Three  Levels  of  Roughness 
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a  power  function  with  a  negative  exponent  of  two.  The  coefficients  of 
these  power  functions  are  related  to  the  RMS  by  means  of  another  power 
expression : 

RMS  =  CA1//2 

where  A  is  the  coefficient  for  a  particular  spectra  and  C  is  a  con¬ 
stant  whose  magnitude  depends  on  the  frequency  range  considered. 

The  results  of  the  simulations  are  shown  in  Fig.  20  in  the  form  of 
RMS  acceleration  of  the  vehicle's  CG  versus  RMS  roughness.  There  is  a 
distinct  response  pattern  for  each  speed,  but  the  influence  of  the  sys¬ 
tem's  nonlinearities  at  the  higher  speeds  on  the  rougher  profiles  caused 
pronounced  deviations  from  a  well-behaved  pattern.  This  is  typical  of 
frequency-dependent  systems;  and,  because  of  this,  several  more  profiles 
and  speeds  should  be  inserted  to  adequately  define  the  response  pattern. 
To  get  some  idea  of  the  response  sensitivity,  an  extra  profile  (RMS  = 
5.0)  was  included  for  the  20-mph  speed  and  two  (RMS  =  4.5,  5.2)  for  the 
10-mph  speed.  The  most  noticeable  aspect  of  the  additional  profiles  was 
a  pronounced  peak  that  occurred  in  the  response  pattern  at  20  mph  when 
the  profile  roughness  changed  from  RMS  =  4.9  to  RMS  =  5.0.  This  is  un¬ 
doubtedly  due  to  the  occurrence  of  resonant  conditions  that  caused  ex¬ 
cessive  contact  between  the  axles  and  the  bump  stops  and  produced  accel¬ 
erations  that  would  be  harmful  to  the  vehicle  as  well  as  its  occupants 
in  an  actual  situation.  Perhaps  for  a  more  realistic  range  of  speeds 
and  terrain  roughness  the  response-roughness  patterns  would  all  be  well 
behaved . 

Since  speed  and  ride  severity  are  of  chief  interest  in  cross¬ 
country  ride  dynamics,  a  more  usable  graph  was  obtained  by  cross 
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Figure  20.  RMS  Acceleration  Versus  RMS  Terrain  Roughness; 
Vertical  Motion  at  Vehicle's  CG; 

M37  Truck  Model 
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plotting  these  results  and  relating  the  vehicle  response  to  vehicle 
speed.  These  results  are  shown  in  Fig.  21.  More  tests  are  needed  at 
speeds  below  5  mph  to  adequately  define  the  initial  portions  of  the 
response-speed  curves  in  Fig.  21.  This  range,  shown  by  dashed  lines, 
would  most  likely  be  the  critical  range  for  an  actual  TVS  system,  partic¬ 
ularly  for  the  rougher  profiles.  The  dashed  lines  were  made  to  inter¬ 
cept  at  a  small  common  value  at  zero  speed.  This  represents  the  small 
RMS  level  due  to  engine  vibration  that  occurs  in  an  actual  vehicle  when 
the  vehicle  is  idle.  These  two  graphs  almost  completely  define  the  TVS 
characteristics  for  this  particular  vehicle  model.  For  example,  for  a 
given  level  of  ride  severity  based  on  RMS  acceleration  and  a  given  level 
of  terrain  roughness  based  on  RMS  elevations,  estimates  of  the  average 
speed  that  the  vehicle  can  attain  and  stay  below  the  selected  ride 
level  can  be  made. 

However,  because  of  the  statistical  nature  of  the  development,  a 
complete  definition  of  the  TVS  characteristics  should  provide  a  deter¬ 
mination  of  the  probability  of  exceeding  the  criteria  of  interest. 

This  is  attained  by  means  of  the  amplitude  probability  distribution 
(APD) ,  which  is  the  probability  that  the  amplitudes  composing  a  wave¬ 
form  will  exceed  a  given  level.  For  digital  data,  an  estimate  of  the 
APD  is  determined  by  simply  counting  the  number  of  points  above  the 
given  level  and  calculating  the  ratio  of  this  number  to  the  total 
number  of  available  points. 

For  the  terrain  profile,  the  APD  serves  as  a  means  for  estimating 
the  occurrence  of  elevations  above  a  given  level.  For  vehicle  response, 
the  APD  serves  to  estimate  the  probability  of  exceeding  a  given  g-level. 
Such  a  distribution  function  was  constructed  from  the  absolute  values 
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Figure  21.  RMS  Acceleration  Versus  Vehicle  Speed; 
Vertical  Motion  at  Vehicle's  CG; 

M37  Truck  Model 
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of  the  CG  accelerations  for  the  test  at  7.5  mph  over  the  terrain  with 
an  RMS  roughness  of  five,  the  results  are  shown  in  Fig.  22.  The  maximum 
acceleration  that  occurred  during  this  test  was  3.85  g.  According  to 
amplitude  distribution,  the  probability  of  exceeding  2.0  g  is  about  0.90, 
while  the  probability  of  exceeding  3.0  g  is  only  about  0.26.  Looked  at 
in  another  way,  when  a  vehicle  is  traversing  a  terrain  of  this  nature  at 
7.5  mph,  it  can  be  expected  to  exceed  3  g  for  about  26  percent  of  the 
duration  of  the  run.  The  response  spectra  will  give  the  frequency  of 
occurrence  of  these  g  levels.  This  leads  quite  naturally  into  a  possi¬ 
bility  for  quantifying  the  intensity  of  a  vibrational  environment  in 
terms  of  human  tolerance.  For  example,  such  quantities  as  the  frequency 
of  occurrence  of  accelerations  above  a  certain  g  level,  or  the  area 
under  the  RMS  acceleration-time  history  curve  (which  is  a  sort  of  im¬ 
pulse  measure  or  measure  of  vibrational  energy)  should  correlate  with 
appropriate  measures  of  human  tolerance. 

Figs.  20,  21,  and  22  thus  completely  define  the  ride  dynamics 
characteristics  for  the  vehicle  model  operating  on  such  random  undu¬ 
lating  profiles.  However,  since  these  profiles  were  constructed  from 
a  basic  spectral  shape  that  has  been  shown  to  represent  most  natural 
and  man-made  surfaces,  it  is  conceivable  that  such  TVS  relations  could 
be  developed  for  any  desired  vehicle  and  serve  to  describe  its  ride 
dynamics  characteristics  for  both  on-road  and  off-road  environments. 

This  graphical  scheme  will  allow  rating  of  vehicles,  on  a  standard 
competitive  basis  of  ride,  as  to  their  capability  to  negotiate  cross¬ 
country  terrains  at  various  speeds.  Such  a  comparison  could  take  the 
form  of  the  hypothetical  relation  in  Fig.  23,  which  shows  the 
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CG  Acceleration,  g 

e  22.  Amplitude  Probability  Distribution  Function  for  CG  Accelerations 
Incurred  at  7.5  mph  Over  Terrain  Profile  with  RMS  Elevation  =  5 


Vehicle  Speed,  mph 


Figure  23.  Comparison  of  Dynamic  Response  of  Three  Vehicles 
Operating  in  a  Giyen  Terrain  Environment 
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response-speed  relations  for  three  vehicles  operating  in  a  given  off¬ 
road  environment.  As  before,  the  complete  evaluation  as  to  which 
vehicle  was  best  suited  to  the  environment  would  require  the  APD  for 
each  vehicle. 

Also,  the  capability  for  constructing  artificial  profiles  to  con¬ 
form  to  any  desired  power  spectrum  will  permit  vehicle  response  evalua¬ 
tions  for  remote  areas  without  physical  testing  of  the  vehicle  in  that 
area.  This  analysis  capability  greatly  enhances  the  utility  of  a 
comprehensive  analytical  model  for  describing  the  overall  mobility 
aspects  of  off-road  vehicles.  It  also  can  serve  as  an  aid  in  the  design 
of  off-road  vehicles  by  permitting  the  study  of  simulated  motions  of 
candidate  vehicles  operating  on  generalized  terrain  types.  The  effects 
of  suspension  components,  for  example,  can  be  studied  in  terms  of  their 
vibrational  isolation  capabilities  by  examining  the  ratios  of  the  RMS 
accelerations  of  the  axles  to  those  at  various  parts  of  the  vehicle’s 
main  frame,  or  in  regard  to  various  independent  suspension  types  by 
studying  the  ratios  of  the  RMS  fore  and  aft  accelerations. 

Application  to  Real-World  Situations 

This  has  been  a  computer  development  of  the  terrain-vehicle-speed 
relations  for  a  computer  model  that  is  representative  of  an  actual 
vehicle  in  the  military's  Table  of  Organization  and  Equipment  (TO&E) 
traversing  computer-constructed  terrain  profiles  whose  power  spectra 
were  representative  of  those  of  actual  profiles.  However,  the  appli¬ 
cation  of  this  scheme  of  analysis  to  actual  vehicles  traversing  actual 
terrains  is  necessary  to  establish  its  success  and  practicality.  The 
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main  problems  that  might  occur  will  probably  be  associated  with  the 
terrain  input.  The  terrain  profiles  in  this  study  were  constructed  to 
fulfill  certain  essential  requirements  of  randomness  and  stationarity . 
Tests  need  to  be  conducted  to  establish  the  degree  of  stationarity 
exhibited  by  representative  terrains.  A  permissible  degree  of  terrain 
stationarity,  which  has  little  effect  on  the  TVS  relation,  needs  to  be 
established.  This  may  turn  out  to  be  a  uniform  standard  applicable  to 
all  vehicles,  or  may  need  to  be  specified  according  to  vehicle  classes, 
such  as  wheeled,  tracked,  etc. 

The  effect  on  both  the  input  and  output  statistics,  as  well  as 
the  frequency  distribution,  of  detrending  (removing  the  longer  terrain 
wavelengths  that  do  not  affect  ride)  the  original  profile  to  achieve 
this  adequate  degree  of  stationarity  needs  to  be  thoroughly  studied. 
These  results  should  establish  the  requirements  necessary  to  develop 
actual  TVS  relations. 
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CHAPTER  IV:  CONCLUSIONS  AND  RECOMMENDATIONS 

Conclusions 

On  the  basis  of  the  results  of  this  study,  it  is  concluded  that: 

a.  Vehicle  response  statistics  obtained  from  traversing 
terrains  with  constant  statistics  and  a  given  terrain 
power  spectra  will  be  similar. 

b .  Three  basic  graphs  can  be  constructed  that  will  completely 
define  the  ride  dynamics  characteristics  for  a  vehicle 
operating  in  both  on-road  and  off-road  environments. 

c ■  This  graphic  scheme  will  rate  vehicles,  on  a  standard 
competitive  basis  of  ride,  as  to  their  capability  to 
negotiate  cross-country  terrains  at  various  speeds. 

d .  The  capability  for  constructing  artificial  profiles  to 
conform  to  any  desired  power  spectrum  will  (1)  permit 
vehicle  response  evaluations  for  remote  areas  through 
simulation  methods,  (2)  enhance  the  utility  of  a  compre¬ 
hensive  analytical  model  for  predicting  overall  mobility 
aspects  of  vehicles,  and  (3)  aid  in  the  design  of  off¬ 
road  vehicles. 


Recommendations 


It  is  recommended  that: 

a.  The  results  of  this  study  be  applied  to  describe  actual 
TVS  relations. 


b.  Tests  with  actual  vehicles  in  real  terrain  environments 
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be  conducted  to  determine: 

(1)  The  degree  of  stationarity  exhibited  by  repre¬ 
sentative  terrains. 

(2)  The  effect  of  detrending  the  original  terrain  to 
achieve  adequate  stationarity. 

(3)  Minimum  length  of  terrain  required  to  insure  adequate 
stationarity  of  the  output  response. 

(4)  The  applicability  of  such  TVS  evaluations  to  compre¬ 
hensive  analytical  cross-country  mobility  models. 
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Abstract 

The  specific  objective  of  this  study  was  to  define  and  formulate 
the  procedures  for  implementing  a  practical,  compact  description  of 
terrain-vehicle-speed  (TVS)  systems  and  thus  establish  a  methodology 
for  relating  measures  of  ride  to  measures  of  terrain  roughness  and 
cross-country  speed.  This  methodology  should  enable  the  determination 
of  reliable  estimates  of  the  average  speed  a  given  vehicle  can  attain 
under  given  vibrational  constraints  as  well  as  the  probabilities  of 
exceeding  certain  specified  levels. 

This  study  was  entirely  a  computerized  effort.  Digital  computer 
simulations  were  conducted  in  which  a  vehicle  was  run  at  several 
selected  speeds  over  terrain  profiles  with  various  levels  of  roughness. 
The  terrain  profiles  were  generated  by  computer  programs  that  con¬ 
structed  and  shaped  sequences  of  random  normal  numbers  to  provide 
representative  profiles  with  a  desired  power  spectrum,  variance,  mean 
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value,  standard  deviation,  and  RMS.  These  programs  also  performed 
the  necessary  statistical  checks  for  normality,  stationarity ,  and 
randomness. 

A  comprehensive,  nonlinear  mathematical  model  of  an  M37  truck 
served  as  the  vehicle.  This  model  allowed  for  the  nonlinearities  in¬ 
herent  in  large  rotational  motions,  the  suspension  elements,  bump 
stops,  loss  of  ground  contact,  and  the  tire  compliance,  which  was 
represented  by  a  cluster  of  radially  projecting  springs. 

The  results  are  presented  in  the  form  of  three  basic  graphs 
comparing  input  and  output  statistics  and  distribution  functions. 

The  input  statistics  consisted  of  terrain  roughness  as  measured  by 
RMS  elevation.  The  output  is  represented  as  the  RMS  vertical  accel¬ 
eration  of  the  vehicle's  center  of  gravity.  A  family  of  curves  was 
developed  relating  vehicle  response  to  terrain  roughness  for  each 
speed  of  vehicle  traversal.  Cross  plots  of  these  established  relations 
provided  a  more  useful  relation  between  vehicle  response  and  speed 
for  various  degrees  of  terrain  roughness. 

The  use  of  such  relations  to  catalog  vehicles  in  terms  of  their 
ride  dynamics  capabilities  is  suggested.  The  TVS  scheme  is  readily 
adaptable  to  computer  models  for  evaluating  and  rating  the  cross¬ 
country  mobility  characteristics  of  vehicles  traversing  large,  non- 
homogeneous  terrain  areas. 
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CHAPTER  VI:  APPENDIX 

APPENDIX  A:  COMPUTER  PROGRAMS 

The  study  discussed  in  this  report  consisted  entirely  of  a  computer 
effort.  The  digital  computer  programs  employed  in  this  study  are  listed 
on  the  following  pages  in  the  proper  sequence  for  developing  the  rela¬ 
tions  established  in  this  thesis.  The  programs  are  written  in  FORTRAN 
IV  programming  language  for  use  on  a  GE-430  computer  system  in  the  time 
sharing  mode.  A  brief  description  of  each  program  is  given  in  the 
following  paragraphs. 

NOISE 

The  first  program  listed,  NOISE,  was  developed  specifically  to 
compute  the  random  terrain  profiles  with  certain  specified  statistics. 

It  first  calls  for  a  name  for  the  output  file  of  terrain  profile  points. 
The  input  variables  are  TAU,  RMS,  ALPHA,  FACT1,  IX,  and  III.  The 
variables  TAU  and  ALPHA  control  the  PSD  shaping  and  it  was  found  that 
setting  the  ALPHA-TAU  product  =  0.055  gave  the  best  normalities  con¬ 
dition  based  on  results  of  the  Chi-square  test.  The  variable  RMS  is 
used  to  determine  the  standard  deviation  of  the  sequence  of  profile 
points  (see  line  No.  1180  in  program  listing).  FACT1  is  used  to  control 
the  RMS  roughness  level  of  the  resulting  profile  points.  For  example,  a 
value  of  0.1168  for  FACT1  produced  a  sequence  with  an  RMS  =  1.0.  To 
obtain  a  profile  sequence  with  an  RMS  roughness  of  4.0  requires  a  value 
for  FACT1  which  is  (4) (0.1168).  The  input  variable  IX  is  the  arbitrary 
starting  point  in  the  random  number  chain;  III  is  the  number  of  desired 
profile  points. 

The  NOISE  program  utilizes  two  library  subroutines:  GAUSS  and 
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RANDU.  The  subroutine  RANDU  generates  a  sequence  of  uniform,  random 
numbers  by  the  conventional  congruence  method.  The  subroutine  GAUSS 
generates  a  sequence  of  normal,  random  numbers  from  the  uniform  se¬ 
quence  by  a  technique  referred  to  as  the  "sum  of  uniform  deviates." 

A  "shift"  routine  adjusts  the  numbers  to  provide  an  absolute  zero 
mean.  The  normal,  random  sequence  is  passed  through  a  numerical, 
low-pass  filter  to  yield  the  final,  proper  frequency  weighted  sequence. 
STANOR 

The  next  program,  STANOR,  provides  all  the  necessary  tests  of 
the  fundamental  assumptions.  The  input  to  this  program  is  the  output 
file  of  profile  points  generated  by  the  NOISE  program. 

This  program  first  performs  the  trend  test  and  run  test  to  check 
the  stationarity  of  the  points  at  the  0.05  level  of  significance.  The 
sequence  of  numbers  is  divided  into  N  equal  intervals  for  the  two 
tests  and  the  mean  value  computed  for  each  interval.  The  resulting 
sequence  of  the  means  is  checked  for  the  appropriate  number  of  zero 
crossings  (run  test)  and  reverse  arrangements  (trend  test).  The  Chi- 
square  test  is  then  conducted  to  check  the  normality  condition.  This 
is  also  accomplished  at  the  0.05  level  of  significance.  For  this  test, 
the  profile  sequence  is  again  divided  into  an  appropriate  number  of 
groups  determined  by  a  formula  that  depends  on  the  total  number  of 
points  in  the  sequence.  Three  statistical  tables,  TABK1,  TABK2,  and 
TABK3,  are  used  in  conjunction  with  this  program. 

PSD 

The  program,  PSD,  first  computes  the  autocorrelation  function  from 
which  it  obtains  the  PSD  estimate  via  the  Fourier  transform  of  the 
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autocorrelation  function.  The  input  file  is  the  file  of  numbers  of 
which  the  PSD  is  to  be  obtained.  Output  files  of  the  correlation 
coefficient  versus  lag  and  PSD  estimates  versus  frequency  are  generated. 
The  basic  input  variables  are:  KKK,  LAG,  DELTAT,  and  FACT,  where  KKK 
is  the  number  of  points  in  the  input  file,  LAG  is  the  maximum  lag  (199 
for  this  study)  used  in  computing  the  autocorrelation  function.  DELTAT 
is  a  value  that  determines  the  spacing  between  the  input  points  and  may 
represent  either  time  or  distance.  FACT  is  a  coefficient  that  multi¬ 
plies  each  input  number  (see  line  No.  1064  in  listing).  This  serves  as 
a  means  of  altering  the  input  numbers  by  a  factor.  If  no  altering  is 
desired,  set  the  value  of  FACT  to  unity. 

The  frequency  w  is  determined  by  the  formula 


hr 

mAt 


h  =  0,1,2, .  .  .m 


where 

m  =  maximum  lags  used  not  to  exceed  199 
At  =  the  equal  intervals  between  points  (may  be  time  or  space 
units) 

This  frequency  u>  is  in  units  of  radians/unit  time  or  distance  and 
must  be  divided  by  2r  to  obtain  frequency  in  cycles/unit  time  or 
distance. 

For  example,  calculate  the  maximum  obtainable  spacial  frequency 
for  m  =  199  and  At  =  1  ft.  The  maximum  frequency  is  obtained  when 


h  =  m  =  199,  i.e. 
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or  in  terms  of  cycle/ft 


TT 

At 


3.14 

1 


=  3.14  rad/ft 


fi 

max 


3.14 

2tt 


3.14 

6.28 


=  0.5  cyc/ft 


It  is  worthy  to  mention  that  the  raw  spectral  estimates  are  smoothed 
with  Hamming  coefficients. 

TRUCK 

The  TRUCK  program,  which  represents  a  two-dimensional  mathematical 
model  for  the  M37,  3/4-ton  truck,  is  a  special  purpose  program  written 
specifically  for  this  study.  The  program  is  composed  basically  of  an 
executive  program  and  two  subprograms;  DIFFEQ ,  which  contains  the  Runge- 
Kutta-Gill  solution  to  the  differential  equations  and  ALGEBR,  which 
performs  all  the  required  algebraic  operations. 

The  input  variables  required  by  the  program  are  clearly  defined 
in  the  program  listing.  However,  some  noteworthy  comments  concerning 
the  segmented  wheel  implementation  are  given  to  enhance  an  under¬ 
standing  of  the  method  used  to  simulate  the  tire-terrain  compliance. 

Each  wheel,  which  represented  a  9.00x16,  8-PR  tire  at  45-psi 
inflation  pressure,  was  divided  into  ten  segments;  five  on  each  side 
of  the  tire's  centerline  as  shown  in  Fig.  Al.  The  measured  load- 
deflection  relation  (Fig.  A2)  for  the  9.00x16  tire  at  45-psi  inflation 
pressure  was  such  that  a  centerline  deflection  of  1.5  in.  required  a 
load  of  2860  lb.  At  this  deflection,  four  spring  segments  are  influ¬ 
enced,  two  on  each  side  of  the  centerline  (Fig,  1).  The  spring  constant 
K  can  be  determined  from  the  statics  equation: 
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F  =  y  K  cos  i  A. 

■  ■>  1  x 

1=1 

where 

A±  =  y  -  THRESH(i)  -  Z  ,  Y±  -  THRESH(i)  -  Z  >  0 

=  0  ,  Y.  -  THRESH(i)  ~  Z  <  0 

=  vertical  height  of  terrain  profile  beneath  i^ 
segment 

Z  =  vertical  displacement  of  axle 
THRESH (i)  =  height  from  the  zero  reference  to  the  i^  spring  of 
the  undeflected  wheel  (see  Fig.  1). 

For  this  case  and  due  to  the  symmetry  of  the  segments  about  the  center- 
line  the  equation  reduces  to 


2 

2860  =  2K  l  1.4  cos 
i=l 


11.7° 

2 


+  0.7  cos  11.7°  + 


11.7° 

2 


where  the  effective  radial  deflections  are: 

Ar  =  Ar  =  1.4  in. 

5  6 

A.  =  A.,  =  0.7  in. 

4  7 

Solving  for  K  yields 

K  =  675  lb/in. 

Defining  GAMMA  =  K  cos  (Jk  =  675  cos  (j>  yields  the  following  relations 
for  the  segments  of  the  front  and  back  wheels: 

GAMMA(l)  =  GAMMA(IO)  =  GAMMA (11)  =  GAMMA(20)  =  411.75 
GAMMA (2)  =  GAMMA (9)  =  GAMMA (12)  =  GAMMA (19)  =  506.25 

GAMMA(3)  =  GAMMA (8)  =  GAMMA(13)  =  GAMMA(18)  =  587.25 

GAMMA (4)  =  GAMMA (7)  =  GAMMA(14)  =  GAMMA (17)  =  648.00 
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GAMMA (5)  =  GAMMA(6)  =  GAMMA(15)  =  GAMMA(16)  =  668.25 
A  similar  relation  is  derived  for  the  threshold  heights  of  each  segment 
THRESH(i)  and  these  relations  are  defined  in  the  program. 

A  mean  spacing  of  3.07  in.  was  determined  to  be  adequate  for  por¬ 
traying  the  projected  spacing  of  the  springs.  As  a  result,  all  profile 
points  were  spaced  3.07  in.  apart  and  no  interpolation  scheme  was 
employed  to  estimate  elevation  between  adjacent  points. 

The  truck  was  then  advanced  point  by  point  along  the  profile  and 
at  each  advance  the  resulting  motions  were  computed. 
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NOISE 


900SNDM 

lOOOSLIB, GAUSS , , ,*** 

1 0 1  oc  ******************************************************** 

1020C  *  LOW-PASS  GAUSSIAN  RANDOM  NUMBER  GENERATOR  * 

1030C  ******************************************************** 

1110  DIMENSION  X(1000),Y(1000) 

1 1  30C  ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 

1132  PRINT, "OUTPUT  FILE  NAME" 

1134  READ  108.XNAME 

1136  108  FORMAT (A  6 ) 

1 1 38  PRINT  ,  "TAU  ,RMS  .ALPHA  .FACT  1  .  IX  .  1 1 1  “ 

1140  READ  ,TAU  ,RMS .ALPHA  .FACT  1  , 1 X , 1 1 1 

1180  SIGMAN:RMS*SQRT(1 . -EXP ( -2. *ALPHA *TAU > ) 

1190  PRINT  102.SIGMAN 

1200  102  FORMAT (8HSIGMAN:  E16.6) 

1 2 1 OC  ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 

1220C  GENERATE  THE  RANDOM  NORMAL  NUMBERS  HAVING  STANDARD 

1 23 OC  DEVIATION  SIGMAN. 

1235C  ******************** 

1249C  ******************** 

1250  AM:0. 

1260  DO  200  I  =  !  .Ill 

1262  SIGSQ=SIGMAN*SI GMAN 

1270  CALL  GAUSS  (IX  .SIGSQ  .AM  ,V) 

1280  200  X(I  ):V 

1281  CALL  SHIFT  (X  .III  .XBAR) 

1282  PRINT  210, XBAR 

1283  210  FORMAT (  "XBAR  AFTER  GAUSS=”E 16.6) 

1285  CALL  STDEVCX  .III  , DEV) 

1286  RMSRADrSQRT (DEV ) *3. 141 5926536*2. 

1287  PRINT  220.DEV.RMSRAD 

1289  220  F0RMAT("DEV:"E16.6,"RMSRAD:“E16.6> 

1290  PRINT  106 

1300  106  FORMAT ("RANDOM  NORMAL  NUMBERS  WITH  STAND.  DEV.  SIGMAN") 
1320  103  FORMAT (E 20 . 10) 

1334C  STANDARD  DEVIATION  IS  COMPUTED  TO  SEE  WHETHER 
1336C  SUBROUTINE  GAUSS (I  X  ,SI GMAN  .AM  ,V  )  GENERATES  CORRECT 
1338C  RESULTS. 

1340  CALL  STDEV(X  .III  , DEV) 

1342  PRINT  330, DEV 

1344  330  FORMAT ("COMPUTED  STANDARD  DEVIATION:  "  E16.8) 

1350C  STANDARD  DEV.  BLOCK  IS  TO  BE  INSERTED. 

1370C  GENERATE  THE  LOW-PASS  GAUSSIAN 

1375  CALL  SHIFT  (X  .III  .XBAR  ) 

1379C  **************************** 

1380  AA:EXP (-ALPHA *TAU) 

138 1C  **************************** 

1390  Y  (  1  )  : 0 . 0 


NOISE  CONTINUED 


1400  DO  300  I =2, III 

1410  300  Y(I  )=X(I  )+Y  (I  -  1  )*AA 
1412  DO  401  1  =  1  ,111 

1414  401  Y(I  )=FACT1*Y(I  ) 

1420  WRITE(2,103)(Y(I ), 1  =  1, III  ) 

1430  CALL  CLOSEF (2,XNAME) 

1435  PRINT  107 

1437  107  FORMAT  (“LOW-PASS  GAUSSIAN  IS  GENERATED  AND  STORED") 

1440  END 

1 99  OC  ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 

2000  SUBROUTINE  RANDU (IX  ,IY ,YFL) 

2010  IY=IX*4099 

2020  IF (IY )5 , 6,6 

2030  5  IY=IY +8 388607+1 

2040  6  YFL=IY 

2050  YFL=YFL*1 . 192093E-7 

2060  RETURN 

2070  END 

208  OC  +++++++++++++++++++++++++++++++++++++++++++++++++++++++ 

2090  SUBROUTINE  STDEVCX  , II I  , DEV) 

2095  DIMENSION  X(1 ) 

2100  SUM  =  0.0 

21 10  DO  420  1  =  1  , 1 1 1 

2120  420  SUM  =SUM  +X (I ) 

2130  SUM1=SUM*SUM 

2140  SUM  2=0.0 

2150  DO  440  1  =  1  , III 

2160  440  SUM2=SUM2+Xd  )*X(I  ) 

2165  AI 1 1 =1 1 1 

2170  SD=SQRT ((SUM  2-SUM  1 /AI 1 1  )/(AIII-l.)) 

2180  DEV=SD 

2190  RETURN 

2200  END 

30  0  OC  +++++++++++++++++++++++++++++++++++++++++++++++++++++++ 

3010  SUBROUTINE  SHIFT(X  ,111  ,XBAR > 

3020  DIMENSION  X(l) 

3030  SUM  =  0 . 0 

3040  DO  300  1  =  1  , 1 1 1 

3050  SUM  =SUM  +X (I  ) 

3060  300  CONTINUE 
3065  AIII =111 

3070  XBAR  =  SUM  /AIII 

3080  PRINT  310, XBAR 

3090  310  FORMAT("THE  TIME  SERIES  IS  SHIFTED  BY  XBAR="E16.8) 

3100  DO  320  1=1,111 

3110  Xd  >=X(I  )-XBAR 

3120  320  CONTINUE 
3130  RETURN 

3140  END 
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GAUSS 


1750C 

1762C . 

1774C 

2006C  SUBROUTINE  GAUSS 
2020C 

2032C  PURPOSE 

2044C  COMPUTES  A  NORMALLY  DISTRIBUTED  RANDOM  NUMBER  WITH  A  GIVEN 

2056C  MEAN  AND  STANDARD  DEVIATION 

2070C 

2102C  USAGE 

2114C  CALL  GAUSS (IX ,S  ,AM  ,V ) 

2126C 

2140C  DESCRIPTION  OF  PARAMETERS 

2152C  IX  -IX  MUST  CONTAIN  AN  ODD  INTEGER  NUMBER  WITH  NINE  OR 
2164C  LESS  DIGITS  ON  THE  FIRST  ENTRY  TO  GAUSS.  THEREAFTER 

2176C  IT  WILL  CONTAIN  A  UNIFORMLY  DISTRIBUTED  INTEGER  RANDOM 

221 OC  NUMBER  GENERATED  BY  THE  SUBROUTINE  FOR  USE  ON  THE  NEXT 

2222C  ENTRY  TO  THE  SUBROUTINE. 

2234C  S  -THE  DESIRED  STANDARD  DEVIATION  OF  THE  NORMAL 
224SC  DISTRIBUTION. 

2260C  AM  -THE  DESIRED  MEAN  OF  THE  NORMAL  DISTRIBUTION 

2272C  V  -THE  VALUE  OF  THE  COMPUTED  NORMAL  RANDOM  VARIABLE 

2304C 

23  ISC  REMARKS 

2330C  THIS  SUBROUTINE  USES  RANDU  WHICH  IS  MACHINE  SPECIFIC 
2342C 

2354C  SUBROUTINES  AND  FUNCTION  SUBPROGRAMS  REQUIRED 

2366C  RANDU 

2400C 

2412C  METHOD 

2424C  USES  12  UNIFORM  RANDOM  NUMBERS  TO  COMPUTE  NORMAL  RANDOM 

24 3 SC  NUMBERS  BY  CENTRAL  LIMIT  THEOREM.  THE  RESULT  IS  THEN 

245 OC  ADJUSTED  TO  MATCH  THE  GIVEN  MEAN  AND  STANDARD  DEVIATION. 

2462C  THE  UNIFORM  RANDOM  NUMBERS  COMPUTED  WITHIN  THE  SUBROUTINE 
2474C  ARE  FOUND  BY  THE  POWER  RESIDUE  METHOD. 

2506C 

25  20C . 

2532C 

2544  SUBROUTINEGAUSS (IX ,S  ,AM  ,V) 

25 5 S  ArO.O 

2570  D050I  =  1 ,12 

2602  CALLRANDU (IX  ,IY ,Y  ) 

2614  IX=IY 

2626  5 OA ; A  +Y 

2640  Vr(A-6.0)*S+AM 

2652  RETURN 

2664  END 
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RANDU 


1750C 

1762C . 

1774C 

2006C  SUBROUTINE  RANDU 
2020C 

20 32C  PURPOSE 

2044C  COMPUTES  UNIFORMLY  DISTRIBUTED  RANDOM  REAL  NUMBERS  BETWEEN 

2056C  0  AND  1.0  AND  RANDOM  INTEGERS  BETWEEN  ZERO  AND 

2070C  2**31.  EACH  ENTRY  USES  AS  INPUT  AN  INTEGER  RANDOM  NUMBER 

2102C  AND  PRODUCES  A  NEW  INTEGER  AND  REAL  RANDOM  NUMBER. 

21  14C 

2126C  USAGE 

2140C  CALL  RANDU  (IX  ,IY  ,YFL) 

2152C 

2164C  DESCRIPTION  OF  PARAMETERS 

2176C  IX  -  FOR  THE  FIRST  ENTRY  THIS  MUST  CONTAIN  ANY  ODD  INTEGER 
221 OC  NUMBER  WITH  NINE  OR  LESS  DIGITS.  AFTER  THE  FIRST  ENTRY, 

2222C  IX  SHOULD  BE  THE  PREVIOUS  VALUE  OF  IY  COMPUTED  BY  THIS 

2234C  SUBROUTINE. 

224 SC  IY  -  A  RESULTANT  INTEGER  RANDOM  NUMBER  REQUIRED  FOR  THE  NEXT 
226 OC  ENTRY  TO  THIS  SUBROUTINE.  THE  RANGE  OF  THIS  NUMBER  IS 

2272C  BETWEEN  ZERO  AND  2**31 

2304C  YFL-  THE  RESULTANT  UNIFORMLY  DISTRIBUTED,  FLOATING  POINT, 

2316C  RANDOM  NUMBER  IN  THE  RANGE  0  TO  1.0 

2330C 

2342C  REMARKS 

235  4C  THIS  SUBROUTINE  IS  SPECIFIC  TO  THE  GE-427 
2366C  THE  SUBROUTINE  SHOULD  PRODUCE  (2**225-1  TERMS 
2400C  BEFORE  REPEATING 

2412C 

2424C  SUBROUTINES  AND  FUNCTION  SUBPROGRAMS  REQUIRED 

2436C  NONE 

2450C 

2462C  METHOD 

2474C  POWER  RESIDUE  METHOD  DISCUSSED  IN  IBM  MANUAL  C20-8011, 

2506C  RANDOM  NUMBER  GENERATION  AND  TESTING 

2520C 

2532C . 

25  44C 

2556  SUBROUTI NERANDU  (IX,IY,YFL) 

2570  IY=IX*4099 
2602  IF  <IY  >5,6,6 
2614  5IY=IY+8388607+t 
2626  6YFL=IY 

2640  YFL  =YFL *  1 . 1 92093E-7 
2652  RETURN 
2664  END 
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STANOR 


1  20  $N  DM 

1 30C  ************************************ 

140C  *  STATI ONARI TY  AND  NORMALITY  TEST  * 

1  5  OC  ************************************ 

160  DIMENSIONX(IOOO) 

165  DIMENSIONZ (50) ,T (50) 

170  DIMENSIONTF (50) 

180  DIMENSI0NTABLE1 (31 ,7) .TABLE 2( 1 5 , 7 ) .TABLE3 (2, 7 ) 

185  COMMONMTEST 

186  MTEST=0 

190  1  OOFORMAT (81  10) 

195  109FORMAT (7F9.3) 

200  101FORMAT (8F10.3) 

210  1 02FORMAT  C8E 1 2. 4) 

220  1 03FORMAT ( 4H  IP =  1 5 ,8E 1 2. 4) 

230  104FORMAT (8F10.5) 

24 OC  TABLE1 (I  ,J  )---P£RCENTAGE  POINTS  OF  RUN  DISTRIBUTION 
25 OC  TABLE3G ,J)---GROUP  NUMBER  DECISION  TABLE 

26 OC  TABLE2CI.J) - PERCENTAGE  POINTS  OF  REVERSE  ARRANGEMENT  DISTRIBU® 

270C  READ  THE  TIME  SERIES 

271  PRINT, "FILE  NAME" 

272  READ839  .XDNAME 

273  839F0RMAT (A  6) 

275  PRINT, "NO.  POINTS,  FACTOR" 

280  READ, III 
290  AIII =111 
300  I P  = 1 000 

310  PRINT745.III 

31 1  745F0RMAT (5H  III =  1 10) 

315  CALLOPENFC 1  .XDNAME) 

318  999 FORMAT (E20. 10) 

320  READG  ,999)<XCI  ), 1=1, III  ) 

321C  FACTOR  IS  MULTIPLIED  FOR  THE  DATA  X444 

360C  *********************** 

370C  READ  THE  PERCENTAGE  POINTS  OF  RUN  TEST  TABLE 
380  CALL0PENF(2,"TABKA 1") 

39  0  READ (2,) ((TABLE  1 (K 1 ,K 2) ,K 2= 1 , 7) ,K  1=  l  ,3  l  > 

440C 

45 OC  READ  THE  PERCENTAGE  POINTS  OF  REVERSE  ARRANGE  DISTRIBUTION  TABLE 

460  CALL0PENF(3, "TABKA  2“ ) 

470  READ(3,)((TABLE2(K 1 ,K 2) ,K 2= 1 , 7) ,K 1=1,15) 

520C  ***************** 

530C  READ  THE  GROUP  NUMBER  DECISION  TABLE 

540  CALLOPENF <4,"TABKA3") 

550  READ(4,)((TABLE3(K l,K2),K2=l ,7),X 1=1,2) 

60 OC  NUMBER  =  NUMBER  IN  SUBSET 

610  NUMBER  = 1 

620C  ****************** 

630  600CONTI NUE 

640C  DETERMINE  NUMBER  OF  GROUP  NG 
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STANOR  CONTINUED 


650  I F  <  200-1 1 1  )339,339,340 
660  340ANG=AIII*(16. 0/200.) 

670  I P : 300 1 
680  GOT0390 

690  3 39 IF (III -20 00)3 49, 35 0,350 
700  350ANG:1 .8 7*(AII I -!.)**( 2.0/5. 0) 

710  I P ; 3002 
720  G0T0390 
730  349CONTI NUE 
740  K2=l 
750  360CONTI NUE 

760  I F  C A 1 1 1 -TABLE3 ( 1 ,K2))361  ,362,363 

770  361ANG=((TABLE3(2,K2+1)-TABLE3(2,K2))/CTABLE3(1 ,K2+1 )-TABLE3( 1 ,K2)) 
780&) 

790&*(AIII -TABLE3C 1 ,K2) )+TABLE3(2,X  2) 

800  I P n 3003 
810  G0T0390 

820  362ANG=TABLE3(2,K2) 

830  I P : 300 4 
840  GO TO 390 
850  363K  2=K  2+ 1 
860  G0T0360 
870  390NG;ANG 
880  PRINT801 ,NG 

885  80 l FORMAT ( l 7H  NUMBER  OF  GROUP=I5) 

89  0  CALLGROUP (X ,111  .ZSTEP  ,NG  ,XMI N ,XMAX ,Z ) 

900  CALLA GROUP (X  .III ,T  ,NG  ,XXMI N .XXMAX ,ASTEP) 

910  CALLTRTEST (T ,NG ,XXMI  N , XXMAX ,TABLE2,IPASS ) 

920  CALLRNTEST (T ,NG ,XXMI N , XXMAX .TABLE  1  .JPASS  ) 

930C  *****  TEST  FOR  NORMALITY  ****************** 

940  CALLCHITNCX .III  ,NG ,XMI  N ,XMAX .ZSTEP ,Z .TF.XBAR ,S ,CHI 2.PCHI ) 

950  PRI NT 40 1 ,PCHI 

960  40 1  FORMAT ( 1 9H  CHI-SQUARE  VALUE  =E12.4) 

970  IF(PCHI -0.05)409,410,410 
980  4 1  OPR  I  NT  420 

990  420FORMAT (30H  NORMALITY  TEST  PASSED  ) 

1000  G0T0499 
101  0  409PRI  NT  430 

1020  430FORMAT (30H  NORMALITY  TEST  NOT  PASSED  ) 

1240  499C0NTI NUE 
1250  END 

1260  SUBROUTI NE GROUP (X ,1 1 1  .ZSTEP ,NG  ,XXMI N  .XXMAX  ,Z  ) 

1265  COMMONMTEST 

1270  DIMENSIONX (2000) ,Z (50) 

1280  104F0RMAT(I5,4£12.4) 

1290  ANG=NG 
1300  AI 1 1 =11 1 
1310  XXMAX  =  - 1 O.E 1 0 
1320  XXMI N  =  1 O.E 1 0 
1330  DO 20 01 s  1  .III 
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STANOR  CONTINUED 


1340  TEMP  1  = AM AX  1 (XXMAX ,X  (I  ) ) 

1350  XXMAX=TEMP1 

1360  TEMP 2=AMI N 1 (XXMI N  ,X  (I ) ) 

1370  200XXMIN=TEMP2 
1380  I P  =  4000 

1390  XXMAX  =XXMAX+1 0.E- 1 0 
1400  XXMIN=XXMIN-I0.E-10 
1420  ZSTEP  =  (XXMAX “XXMI N )  /ANG 
1430  TEMP  1 =XXMI N 
1440  TEMP2=TEMP 1+ZSTEP 

1445  IF (MTEST.NE.0)G0T0444 

1446  107F0RMAT (/,6H  GROUP) 

1449  444CONTINUE 

1450  D0300N=1 ,NG 
1460  Z(N)=0.0 
1470  D0310I  =  1 ,1 1 1 

1480  I  F(X  (I  )-TEMP  1)310,310,31  1 

149  0  3 1 1 1 F (TEMP 2-X (I ))310,310,315 

1500  315Z<N)=Z(N)+1.0 

1510  31 OCONTINUE 

1520  I P  =  4008 

1530  AN  =  N 

1535  IF(MTEST.NE.0)GOTO445 
1540  PRINT104.N, TEMPI, TEMP2,Z(N) 

1545  445C0NTI NUE 

1550  TEMPUTEMP2 

1560  300TEMP 2=TEMP 2+ZSTEP 

1565  MTEST =MTEST+1 

1570  RETURN 

1580  END 

1590  SUBROUTI NETRTEST (Z .NG.XXMIN .XXMAX ,TABLE2,IPASS ) 

1600  DIMENSIONZ (50) ,A(50) ,TABLE2( t 5 , 7) 

1610  1 04F0RMAT ( 4H  IP  =1  5 ,8E 1 2. 4) 

1615  108FORMAT(I5,F10.0) 

1620C  *************  ^REVERSE  ARRANGEMENT  COMPUTATION 

1625  PRI NT  1 09 

1626  1 09FORMAT(/,21H  REVERSE  ARRANGEMENT  ) 

1630  NGM 1 rNG-  1 

1640  ANG;NG 

1650  DO  200N= 1 , NGM  1 

1660  A(N)=0.0 

1670  NP1=N+1 

1680  D0210J=NP 1  ,NG 

1690  IF(Z(J)-Z(N)) 220 ,210,210 

1700  220A (N ) =A (N )  + 1 .0 

1710  21 OCONTI NUE 

1730  PRI  NT  1  08  ,N  ,A  (N  ) 

1740  200CONTI NUE 

1750C  *******  TOTAL  OF  REVERSE  ARRANGEMENT 

1760  TREVrO.O 
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STANOR  CONTINUED 


1770  DO 25  ON  :  1  ,NGM  1 
1780  250TREV:TREV+A  (N) 

1790  PRINT773 

1792  773F0RMAT(/, "SUMMARY  OF  TREND  TEST") 

1794  PRI NT  774 ,TR£V 

1796  774F0RMAT ( '  TOTAL  REVERSE  ARRANGEMENT  :  "E12.4) 

1810C  NOTE  THAT  IT  IS  TEMPORARY  ASSUMEDE  THAT  K2L=3,AND  K2U=5  IN  THE 
1820C  TABL£2(I ,J)  ARRAY 

1830C  *  *  ^COMPUTE  THE  LOWER  LIMIT  BY  LINEAR  INTERPOLATION.  (ASSUME  B 

1840C  IS  GREATER  THAN  10  AND  LESS  THAN  100) 

1850  IF(NG-10)309,309,301 
1860  301IF(100-NG)309,309,302 
1870  309PRI NT308 ,NG 

1880  308FORMAT ( 48H  NG  IS  GREATER  THAN  100  OR  LESS  THAN  10.  NG  =110) 
1890  STOP 
1900  302C0NTI NUE 
1910  D0641K 1:2,15 

1915  IF (ANG.GT.TABLE2(K 1 , l))G0T0641 
1920  IFCANG.EQ .TABLE2CK 1,1))G0T0642 
1925  XX 1 :TABLE2(K 1-1,1) 

1930  XX2:TABLE2(K 1,1) 

1935  Y 1=TABLE2(K 1-1 ,3) 

1936  Y2:TABLE2(K  1  ,3) 

1940  ALOW: (Y  2-Y 1 )*(ANG-XX l ) /(XX  2-XX 1 )+Y 1  / 

1945  YY 1 :TABLE2(K 1-1,6) 

1946  YY2:TABLE2(K  1,6) 

195  0  AUP: (YY  2-YY 1 )*(ANG-XX 1 ) / (XX 2-XX 1 )+YY 1 
1953  GOTO 643 

1955  642AL0W:TABLE2(K 1 ,3) 

1960  AUP:TABLE2(K 1 ,6) 

1965  GOTO 643 
1970  64 1C0NTI NUE 
1975  643C0NTI NUE 
2034  PRINT775.AL0W 

2036  775F0RMAT ( "  LOWER  LIMIT  OF  REVERSE  ARRANGEMENT  =  "E12.4) 

2038  PRI NT  776, AUP 

2040  776F0RMAT ( "  UPPER  LIMIT  OF  REVERSE  ARRANGEMENT  :  'E12.4) 

2045  5 1 OFORMAT ( 6H  ALOW:E 1 2. 4 , 5H  AUP=E12.4,6H  TREV:E12.4) 

2050C  TEST  WHETHER  ANG  IS  BETWEEN  THE  LIMIT 

2060  IF(TREV- ALOW) 39 0,400, 400 

2070  4001 F (AUP -TREV ) 390 , 41 0 , 4 1 0 

2080  39 01  PASS :0 

2090  PRINT391 

2100  39 1  FORMAT ( "  TREND  TEST  IS  NOT  PASSED") 

21 10  GOTO  420 
2120  4 1 01  PASS : 1 
2130  PRI NT 4 1 1 

2140  4 1  1  FORMAT ( "  TREND  TEST  IS  PASSED  ') 

2150  420CONTI NUE 
2160  RETURN 
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STANOR  CONTINUED 


2170  END 

2180  SUBROUTI NEA GROUP (X  ,1 1 1  ,T ,NG ,XXMI N  .XXMAX , ASTEP  ) 

2I90C  THE  PURPOSE  OF  AGROUP  IS  TO  GROUP  THE  SERIES  X(I)  AND  AVERAGE  I 

2200C  AND  STOR  THEM  IN  T(I)  ARRAY.  NOTICE  THE  DIFFERENCE  BETWEEN  GROB 

221 OC  AND  AGROUP. 

2220  DIMENSIONX (2000) ,T(50) 

2230  1 04FORMAT ( 4H  IP =1 10 ,8E 1 2. 4) 

2240  AII1  =111 
2250  ANGrNG 

2260C  COMPUTE  THE  INTERVAL  IN  FLOATING  MODE 

2270  TEMP  1  =  A 1 1 1 /ANG 
2280  ASTEP  =TEMP  1 

229 OC  THESE  CODES  MAY  BE  TAKEN  OUT  *  *  *  *  * 

235  OC 

2360  1 52F0RMAT ( 7F8. 2) 

2370C  ******************* 

2380  T£MP2=XXMIN 
2390  D0200N: 1 ,NG 
2400  TCN)=0.0 
2410  AN  =N 
2420  NT  =  0 

2430  TEMP  3=TEMP I *(AN- 1 . 0 ) 

2440  TEMP4=TEMP1*AN+1.0E-10 
2450  DO 2 1 01  =  1 .III 
2460  AI:I 

2470  IF(AI-TEMP3)219, 219,220 
2480  2201 F (TEMP 4-AI  )219,219,230 
2490  230T(N)=T(N)+X(I  ) 

2500  NT  =NT  +  1 
2510  219C0NTINUE 
2520  21 OCONTI NUE 

2530C  *******  print  THE  INTERMEDIATE  STEP 

2560C  ******* 

2570  ANT  =  NT 

2580  T (N )=T (N ) /ANT 

2610  200CONTI NUE 

2615  PRINT310 

2616  3 1 OFORMAT ( / , 23H  AVERAGED  TIME  HISTORY  ) 

261 7  D031 IN:  1  ,NG 

2618  31 1PRINT31 2.N ,T (N ) 

2619  312F0RMAT (I 5.E12.4) 

2620  RETURN 
2630  END 

2640  SUBROUTI NERNTEST (X ,NG ,XXMI N .XXMAX .TABLE  1 .JPASS ) 

265 OC  THE  MAIN  PROGRAM  CALL  RNTEST (Z ,NG ,XXMI N .XXMAX .TABLE  1 .JPASS ) 

2660  DIMENSIONX (50), TABLE  1(31 ,7),X 1(50) 

2670  1  04FORMAT ( 4H  IP =1 5 , 4E 1 2. 4 ) 

268 OC  *********COMPUTE  MEAN  VALUE 

2683  189C0NTI NUE 

2684  18  7F0RMATU  5 ,3E  1  2.  4) 
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2690  TEMP  =  0.0 

2700  DOSOON=I ,NG 

2710  800T£MP:TEMP+X(N) 

2720  ANG:NG 

2730  XMEAN:TEMP/ANG 

2760C  ********DETERMINE  NUMBER  OF  RUN 
2770  NRUN=t 

2780  I F (X (  1 ) -XMEAN )  267 , 268 , 268 

2790  2671 CHG=- 1 

2800  G0T0270 

2810  2681 CH  G=  1 

2820  270CONTI NUE 

2830  DO  2601 : 2 ,NG 

2840  IF  (X  <1  ) -XMEAN) 26 2, 261 ,261 

2850  26  IK  =  1 

2860  G0T0263 

2870  262X=-1 

28  8  0  GOTO 263 

2890  263IF(ICHG-K >265,260,265 
2900  265NRUN  =NR  UN  +  1 
2910  ICHG=K 
2920  260CONTI NUE 

2950C  *******COMPUTE  LOWER  LIMIT  BY  LINER  INTERPOLATION 

2960C  (ASSUME  NG/2  IS  GREATER  THAN  5  AND  LESS  THAN  100) 

2970  NNG-NG/2 

2980  IF (NNG-5)309,309,301 

2990  301IF(IOO-NNG)309,309,302 

3000  309PRI NT 308  ,NNG 

3010  308FORMAT ( 45H  NNG  IS  GREATER  THAN  100  OR  LESS  THAN  5.  NG=I10) 

3020  STOP 

3030  302CONTI NUE 

3040  ANNG=NNG 

3050  DO 64 1 K 1:2,31 

305  5  IFCANNG.GT.TABLEKK  1  ,1))G0T0641 
3060  IF (ANNG.EQ. TABLE 1(K 1,1))G0T0642 
3065  XX UTABLEl  (K  1-1 , 1) 

3070  XX2:TABLE1(KI,1) 

3075  Yt=TABLEl(Kl-l,3) 

3080  Y2=TABLE1(K 1,3) 

3085  ALOW: (Y 2-Y 1 )*(ANNG-XX 1 )/(XX2-XX 1 )+Y 1 
3090  YY  UTABLE1  (K  1-1  ,6) 

3095  YY  2=TABLE 1 (K 1 ,6) 

31 00  AUP  =  (YY  2-YY 1  )*(ANNG-XX 1 )/(XX2-XX 1 )+YY  1 
3105  GOTO 643 

3110  642AL0W=TABL£1(K 1 ,3) 

31 15  AUP  =TABLE 1 (K 1 ,6) 

3120  G0T0643 
3125  641CONTINUE 
3130  643C0NTI NUE 
3174  PRI NT  790 
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3176  790F0RMATC/, "SUMMARY  OF  RUN  TEST") 

3178  PRINT791  ,NRUN 

3180  79 1  FORMAT ( "  TOTAL  NUMBER  OF  RUN  =  "I10> 

3181  PRINT793.AUP 

3182  PRINT  792, ALOW 

3184  792F0RMAT  C "  LOWER  LIMIT  OF  RUN  ="£12.4) 

3185  5 1 OFORMAT ( 6H  ALOW=E l 2. 4 ,5H  AUP=E12.4,6H  NRUN=E12.4) 

3186  793F0RMAT ( "  UPPER  LIMIT  OF  RUN  =  “E12.4) 

3190C  TEST  WHETHER  ANNG  IS  BETWEEN  THE  LIMITS 
3200  IF (ANNG- ALOW >390,400, 400 

3210  40 01 F(AUP -ANNG) 39 0,410, 410 
3220  390IPASS=0 
3230  PRI NT 39 1 

3240  39 1  FORMAT  < “  RUN  TEST  IS  NOT  PASSED  ") 

3250  GOTO 420 
3260  4 1 01  PASS = 1 
3270  PRI NT 4 1  1 

3280  4 1 1  FORMAT  < "  RUN  TEST  IS  PASSED  “) 

3290  420CONTI NUE 
3300  RETURN 
3310  END 

3320  SUBROUTINECHITN(X  ,N  ,NG  ,XMI N ,XMAX .STEP ,F  ,TF ,XBAR  ,S  ,CHI  2.PCHI  ) 
3330  DIMENSIONF ( 1 ) ,TF ( 1 ) ,X ( 1 ) 

3340  2F0RMAT (4X 1 0FI0.6) 

3350  SX  =  0. 

3352  PRI NT  79 1 

3354  791F0RMATC/, "SUMMARY  OF  NORMALITY  TEST") 

3360  SX2=0. 

3370  DO 2001  =  1  ,N 
3380  SX=SX+X<I ) 

3390  SX2=SX2+X  (I  >*X  (I  ) 

3400  200CONTI  NUE 
3410  XBAR=SX/N 

3420  S=SQRT ( (SX 2- (SX *SX /N ) )/(N- 1 ) ) 

3430  CALLGROUP  (X  ,N  .STEP  ,NG  ,XMI N  ,XMAX  ,F) 

3440  X 1  - (XMI N -XBAR ) /S 
3450  TEMP =  STEP /S 
3460  D021  01  =  1  ,NG 
3470  X2=Xl+TEMP 
3480  A  1  =RNORM (X  1  ) 

3490  A 2=RN0RM (X 2 ) 

3500  TFCI )  =  (A  2-A 1 )*N 

3510  PRINT2.X  1  ,X2,A  1  ,A2,TF(I  ) 

3520  X  1 =X 2 
3530  21 OCONTI NUE 
3540  CHI  2=0. 

3550  DO 2201  =  1  ,NG 

3560  CHI 2=CHI 2+( CF (I )-TF(I ))**2)/TF(I ) 

3570  220C0NTI NUE 
3580  DF  =NG-  1 
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3590  PCHI=CHIPR(DF,CHI2) 

3600  RETURN 
3610  END 

3620  FUN C TI ONR NORM  (Z  ) 

3630  DIMENSIONA (6) 

3640  A(l)=4.9867347E-2 

3650  A(2)=2.1 1410061E-2 

3660  A <3)=3.2776263E-3 

3670  A (4)=3.80036E-5 

3680  A  (5)=4.88906E-5 

3690  A(6)=5.383E-6 

3700  I F (Z )100,l 10,120 

3710  100C0NTINUE 

3720  ASSIGN 160T0JEL 

3730  Z 1=-Z 

3740  GOTO  130 

3750  110C0NTINUE 

3760  B=  .5 

3770  GOTO  1 60 

3780  1 20CONTI NUE 

3790  ASSIGN 150TOJEL 

3800  ZUZ 

3810  1 30C0NTI NUE 

3820  B:1 . 

3830  DO  1 401  =  1 , 6 
3840  B  =  B+A  (I  )*Z  1**1 
3850  1 40C0NTI NUE 
3860  Bs.5/(B**16) 

3870  GOTOJEL , ( 150, 160) 

3880  1 50C0NTI NUE 
3890  B= 1  .-B 
3900  1 60CONTI NUE 
3910  RNORM=B 
3920  RETURN 
3930  END 

3940CCHIPR  CHISQUARE  PROBABILITY 
3950  FUNCTIONCHIPR CDF.CHSQ  ) 

3960  A=.5*DF 
3970  X=.5*CHSQ 
3980  IF (X >90,90, 100 
3990  90CHIPR  =  1  . 

4000  GOTO  1 70 
4010  100TERM:!. 

4020  SUM : 0 . 

4030  COFN=A 

4040  I F ( 1 3 . —X  >110,110,1 20 
4050  1 1 01 F(A-X ) 1 40 , 1 40 , 1 20 

4060C  CONVERGENT  SERIES  FOR  X  .LT.  A  OR  .LT. 
4070  120C0N=1. 

4080  FACT :-A 
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4090  1 30TEMP  =SUM 

4100  SUM:SUM+TERM 

41 1 0  COFN=COFN+l 

4120  TERM=TERM*X/COFN 

4130  I F (SUM -TEMP ) 1 60 , 150,130 

4140C  ASYMPTOTIC  SERIES  FOR  X  .GTE.  A  AND  X  .GTE.  13. 

4150  1 40CON -0 , 

4160  FACT  =X 

4170  1 50TEMP  =SUM 

4180  SUM  =SUM+T ERM 

4190  C0FN=C0FN-1 

4200  RATI 0=C0FN /X 

4210  TERM=TERM*RATIO 

4220  I F (SUM -TEMP  >160,160,150 

4230  1 60CHI PR  =C0N+EXP (ALOGCSUM )-X+A*ALOG(X )-ALOGC GAMMA (A ,0. ) ) )/FACT 
4240  1 70C0NTI NUE 

424 1C  ***************************************** 

4247C  **************************************** 

4248  RETURN 
4250  END 

5000  FUNCTI ONGAMMA (U ,V) 

5010  DIMENSIONS (8) 

5020  B(1 )=. 035868343 

5030  B<2)=-. 193527818 

5040  B<3)=. 482199394 

5050  B  (4):-. 756704078 

5060  B<5)=. 918206857 

5070  B(6)=-. 897056937 

5080  B(7)=. 988205891 

5090  B(8):-. 577191652 

5100  A=U 

51  10  X=V 

5120  XN=15.0 

5130  1 I F (A ) 4 , 2, 3 

5140  2IFCX-1. 0)4000, 4000, 5000 

5150  4Z  =AI NT (A  ) 

5160  Y=ABS(Z-A) 

5170  EPS=3.0E-8 

5180  I F (Y-EPS ) 21 ,21 ,5 

5190  21A=Z 

5200  GOTO 22 

5210  5Y  = 1 .0-Y 

5220  IF (Y-EPS )6,6,3 

5230  6A=Z-1 .0 

5240  22IFCX-1. 0)4000, 4000, 5000 

5250  3IFCX )7, 1000,7 

5260  7IF(X-1.0)9,9,8 

5270  8IFCA )5000, 10, 10 

5280  1 01 F((X /A )-l. 0)1000, 5000, 5000 

5290  9IF(ABS(A)-l.)ll, 1000, 1000 


STANOR  CONTINUED 


5300  I1IFCX-. 1)1000, 12, 12 

5310  1 21 FCX-. 2)13, 14,14 

5320  13XN=140.0 

5330  GOT05000 

5340  14IF(X-.4)15,16,16 

5350  15XN=80.0 

5360  GOT05000 

5370  1 6IF (X- .6) 1 7, 18 , 18 

5380  17XN=60.0 

5390  GOT05000 

5400  181 F(X- .8 ) 19 , 20 , 20 

5410  19XN=40.0 

5420  GOT05000 

5430  20XN=20.0 

5440  5000XFT  =X 

5450  N=XN 

5460  DO 50011:1  , N 

5470  XFT  =  CXN  AFT)  +  1  .0 

5480  XFT=((XN-A)/XFT)+X 

5490  5001XN=XN-l .0 

5500  TEM  = CALOG (X ) *A )-X 

5510  X FT: EXP CTEM  ) /XFT 

5520  ANS :XFT 

5530  GOT06000 

5540  4000SUM =  0.0 

5550  EPS : 1 4 OE -8 

5560  XM :  1  .  0 

5570  XMT:1 .0 

5580  EX  =  -  1 . 0 

5590  TEM=X 

5600  4001Y=TEM/(XMT*XM) 

5610  401  1SUM=SUM+EX*Y 

5620  IF  CABS (Y ) -EPS ) 4003 , 4002 , 4002 

5630  4002TEM=TEM*X 

5640  XM=XM+1.0 

5650  XMT=XMT  *XM 

5660  EX=-EX 

5670  GOTO 400 1 

5680  4003E  =  -ALOGCX)-. 5772 15 66-SUM 
5690  IFCA )4005, 4004, 4005 
5700  4004ANS :E 
5710  GOTO 6000 

5720  4005IFCA+2. 0)4010, 4009, 4008 
5730  4008ANS:-E+EXP(-X )/X 
5740  GOT06000 

5750  4009SUM:( 1 ,0/X )-( 1 .  0/(X*X  )  ) 
5760  ANS=(E-EXP C-X)*SUM)*.5 
5770  GOT06000 
5780  40 1 OEX  =  -1.0 
5790  N=-A-1.0 
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5800  XMT:1 .0 
5810  SUM: 0 . 0 
5820  TEM :X 
5830  D040061 = i ,N 
5840  TEM:TEM*X 
5850  XM:I+1 
5860  Y=XMT/TEM 
58  70  SUM:SUM+EX*Y 
5880  EX  :  -  EX 
5890  4006XMT  =XMT *XM 
5900  SUM :SUM+1 . 0/X 
5910  Z  -  E  -EXP  (  -X  )  *S  UM 
5920  Y :ABS  (A  ) 

5930  EX  :  - 1 . 0 
5940  XMT  =  1  .0 
5950  N=Y 

5960  DO  40071  : 2 ,N 
5970  XM:I 
5980  XMT  :XMT  *XM 
5990  4007EX:-EX 
6000  ANS  :Z  *EX /XMT 
6010  GOTOSOOO 

6020  1000IFCA-1 .0)  1001 , 1003, 1004 

6030  1004IFCA-2. 0)1003, 1003, 1002 

6040  1  003AT :A-  1 . 0 

6050  ATM :B ( 1 ) *AT 

6060  DO  10051:2,8 

6070  1 005ATM: (ATM+8 (I ))*AT 

6080  ANS:ATM+1.0 

6090  G0T02000 

6100  1001IFCA >1007, 1006, 1006 

61  1  0  1  006BN : A 

6120  BS : A 

6130  GOTO  1010 

61  40  1 0  0  7B  N : A 

6150  BS:A 

6160  101 1IF(BS)1  009, 1010, 1010 

61  70  1 009BS :BS  +  1 . 0 

6180  BN:BS*BN 

6190  GOTO  1011 

6200  1 0 1 OATM :BS  *B ( 1 ) 

6210  DO  1 01 21 =  2,8 

6220  1 0 1 2ATM: (ATM+B (I ))*BS 

6230  ANS:(ATM+1 .0)/BN 

6240  GOT02000 

6250  1 002BN : A- 1,0 

6260  BS :BN 

6270  1013IF(2.0-BS)I015, 1014, 1014 
6280  1 0 1 5B5 :BS -1.0 
6290  BN:B5*BN 
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6300  GOTO  1013 
6310  I014BS=BS-1 .0 
6320  ATM=BS*B(1) 

6330  DO  1 01 61  =  2,8 

6340  1  0 1  6ATM  :  (ATM+B  (I  ))*BS 

6350  ANS=(ATM+1 .0>*BN 

6360  2000IF(X)2001 ,6000,2001 

6370  200 1W= 1 . 0 

6380  XM=1 .0 

6390  SUMrl .0 

6400  EPS: 1 , 0E-8 

6410  2002W: (X /(A+XM ) )  *W 

6420  SUM=SUM+W 

6430  XM:XM+1.0 

6440  I F (ABS (W ) -EPS ) 2003 , 2002 , 2002 
6450  2003T=A*ALOG(X) 

6460  T=EXP (T-X)/A 
6470  ANS:ANS-(T*SUM) 

6480  6000GAMMA :ANS 
6490  RETURN 
6500  END 
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900SNDM 

1 OOOC  ***************************************************** 

1010C  *  AUTO-CORRELATION  USING  FILED  DATA  * 

1020C  ***************************************************** 

1030  DIMENSION  X(IOOO) 

1031  DIMENSION  RX (200) ,PX (200) ,SPX (200) 

1032  PRINT  , 'INPUT  FILE,  AUTO-CORR.  PLOT  FILE,  PSD  PLOT  FILE" 

1034  READ  910,  XNAME 1 .XNAME2.XNAME3 

1036  910  FORMAT (3A  6) 

1038  PRINT, "KKK, LAG, DELTAT, FACT" 

1040  READ.KKK  .LAG, DELTAT, FACT 

1042  LLAG=LAG+1 

1050  CALL  0PENF(1, XNAME 1) 

1060  READ(  1  ,920)(X(K),K  =  1,KKK) 

1061  920  FORMAT (E20. 10) 

1062  DO  144  K  s  1  ,KKK 
1064  144  X(K)=X(K)*FACT 
1070  PRINT  110 

1080  110  FORMAT (20H  RAW  TIME  SERIES  ) 

1095  104  FORMAT ( 5E 1 2, 4) 

1 1 OOC  COMPUTE  THE  MEAN  VALUE  AND  SHIFT  THE  DATA  TO  MEAN  ZERO 
1110  SUMsO.O 
1120  DO  150  K  s  1  ,KKK 
11  30  150  SUM=SUM+X(K) 

1140  AKKK  =KKK 
1150  XBAR=SUM/AKKK 

1155  PRINT  107.XBAR 

1156  107  FORMAT (6H  XBAR=E12.4) 

1160  DO  160  K  =  1  ,KKK 

1170  160  X(K)=X(K)-XBAR 
1180  PRINT  103 

1190  103  FORMAT (32H  TRANSFORMED  INPUT  DATA  ) 

1210  NPOINT  =KKK 

1220  DO  300  I  si  ,LLAG 

1230  SX=0.0 

1240  M  =  I-1 

1250  NX sNPOI NT  -M 

1260  FNX  =NX 

1270  DO  305  J  =  1  ,NX 

1280  K=M+J 

1290  305  SX=SX+X(J)*X(K  ) 

1300  300  RX (I )=SX/FNX 
1310  PRINT  400 

1315  400  FORMAT ( 20H  AUTO-CORR.  FT.  ) 

1330  DO  401  Isl,LLAG,4 
1332  J=I 
1334  J3=I+3 

1340C  PRINT  222,((J,RX(J)),J=I  ,J3) 

1349  401  CONTINUE 

1350  222  F0RMAT(4(I 4 ,E 1 2. 4 ) ) 
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1352  PRINT  520 ,RX ( 1 ) 

1354  520  FORMAT ( 7H  RX ( 1 ) =E 1 2. 4) 

1355  XNAME=XNAME2 

1400C  *******poWER  SPECTRAL  DENSITY  COMPUTATION******** 

1405  DELFr 1 . /(DELTA T*AKKK ) 

1410  PI=3. 1415926536 

1420  FLAGrLAG 

1430  Q=PI  /FLAG 

1440  CONST  =2.*DELTAT /PI 

1450  RXCl )=0.5*RX(1 ) 

1460  RX (LLAG) =  0. 5*RX (LLAG) 

1470  Sl=-Q 

1480  DO  1401  I H  =  1  .LLAG 

1490  PX(IH)  =  0.0 

1500  S  1  =S  1-H3 

1510  S=-S1 

1520  DO  1402  JP  =  1  .LLAG 

1530  S=S+S1 

1540  1402  PX(IH)=PX(IH)+RX(JP)*COS(S> 

1550  1  401  PX(IH)=PX(IH)*CONST 

1560  SPX(1 )=0.54*PX(1)+0.46*PX(2> 

1570  SPX (LLAG)=0.54*PX (LLAG)+0.46*PX (LLAG-1 ) 

1580  KK=LLAG-1 
1590  DO  1415  J - 2,KK 

1600  1415  SPX(J)=0.54*PX(J)+0.23*(PX(J+1)+PX(J-1)) 

1602  WRITE(5,920)(SPX(I ),I=1 .LLAG) 

1610  PRINT  500 

1620  500  FORMAT (30H  POWER  SPECTRAL  DEN.  ) 

1625  XNAME=XNAME3 

1626  I  DEV -  5 

1627  CALL  PLOT 4 (SPX .LLAG .XNAME ,1  DEV ) 

1640  DO  550  I  -  l  .LLAG 

1641  AI=I 

1642  TEMP 13:AI *DELF 

1650  PRINT  554,1  .TEMP  13, SPX(I ) 

1676  550  CONTINUE 

1677  554  FORMAT (I 5,E12.4,E16,6) 

1680  PSUMrO.O 

1690  DO  510  1  =  1  .LLAG 

1700  510  PSUM =PSUM+SPX (I )*DELF 
1710  PRINT  511  ,PSUM 

1720  511  FORMAT (6H  PSUM=E12.4) 

1999  END 

2000  SUBROUTINE  PL0T4CX ,KKK .XNAME  ,IDEV) 

2010  DIMENSION  X(l) 

2015  DIMENSION  Y (71 ) ,ISYM(71) 

2020  DATA  IBLANK/3H  / 

2030  DATA  MINUS/3H-  / 

2040  DATA  IAI/3HI  / 

2050  DATA  ID0T/3H.  / 
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PSD  CONTINUED 


2055  558  FORMAT (A  6) 

20S0C  FIND  THE  MAXIMUN  VALUE  IN  ARRAY  X(I) 

2070  TEMP  1  =  0,0 

2080  DO  iOO  K=1,KKK 

209  0  TEMP  2=ABS (X (K ) ) 

2100  YAXISL:  AMAX 1 (TEMP  1 ,TEMP 2) 

21  10  300  TEMPUYAXISL 

2115  YAXISL=YAXISL+YAXISL/35. 

2140  DELY  =  YAXISL/35. 

2150  DO  200  I Y= I ,71 

2160  AIY=IY 

2170  Y(IY)  =  -YAXISL+DELY*(AIY-1.) 

2190  200  CONTINUE 

220 OC  PRINT  215, (Y (IY )  ,1 Y  =  1 ,71 ) 

2210  215  FORMAT (  5E 1 2. 4) 

2220  202  FORMAT ( "***************************************“ 


2230  PRINT  202 

2235  WRI TE (I  DEV , 202) 

225 OC  INITIAL  PLOT  SET 
2260  DO  401  J  =  1 ,71 

2270  401  ISYM(J)=MINUS 
2280  ISYM (35 ) =IAI 

2290  DO  402  K  =  1  ,KKK 

2300  DO  403  J=2,71 

2310  IFCY(J).GT.X(K).AND.Y(J-1).LT.X(K))  GO  TO  410 
2320  GO  TO  403 
2330  410  ISYM  (J- 1 ) =1  DOT 
2335  403  CONTINUE 

2340C  PRINT  420 , (ISYM <1 > ,1  =  1 , 7 1 > 

2350  420  FORMAT (71A1) 

2355C  +++++++++++++++++++++++++++++ 

2356  WRITE  (I  DEV;  420)  (ISYM  (I  )  ,1  =  1, 71) 

235  7C  ++++++++++++++++++++++++++++++ 

2360  DO  440  1=1,71 
2370  440  I SYM  (I  )  =  IBLANK 

2380  ISYM(35)=IAI 

2381  402  CONTINUE 

2382  PRINT  558.XNAME 

2391  PRINT202 

2392  WRI TE (I  DEV ; 202) 

239 3C  ++++++++++++++++++++++++++++++ 

2394  CALL  CL0SEF(5,XNAME) 

239  5C  ++++++++++++++++++++++++++++++ 


2399  RETURN 

2400  END 


) 
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TRUCK 


1$LIB ,DI FFEQ 
2SLIB .ALGEBR 
3$RPC 
4$N  DM 
5STTY  ,120 

100  COMMON  NSTEPS .DELTAT ,A ,B .MASS .MASS  1 .MASS  2.CP0S 1 ,CPOS  2.CNEG 1 .CNEG2 
110  COMMON  C 1 1 ,C21 ,1 NRTI A .FORCW 1 .F0RCW2 .SPDEF1 .SPDEF2 .DSPDF 1 .DSPDF2 
120  COMMON  FORCK 1 .FORCK  2.VAR 1 ,VAR2,VAR3,VAR4,\/AR5,VARS,VAR7,VAR8 

130  COMMON  ZETA  1  .ZETA2.DZETA2.ETA  1  .ETA2.DETA2.NU  1  .NU2.DNU2 
140  COMMON  MU1.MU2.DMU2 

150  COMMON  THRESH (20) .GAMMA (20) .SEGDEF (20) ,Y (47) 

160  REAL  NU  1  .NU2.MU 1 .MU2.MASS .MASS  1 .MASS 2, 1 NRTI A 
165  IBELL-458752 

170  PRINT  1 

180  1  F0RMAT(18X,33(1H*),/,18X,1H*,31X,1H*,/, 

1 9 0 &  18X.33H*  MURPHY  'S  M-37  BACK-UP  MODEL  *,/, 

2004  18X,1H*,31X,1H*,/,18X,33(1H*),///) 

21  OC - DATA  INITIALIZATION - 

220  A=64.8 

230  B: 47. 2 

240  I NRTI A  =  25446. 

250  MASSES. 05 

260  MASS  1  =  .94 

270  MASS  2:. 8 

280  CPOS  1=11.8 

290  CNEG 1=23.8 

300  CPOS  2= 1 0. 2 

310  CNEG2:44. 

320  THRESH <1  )  =  7. 02 

330  THRESH (2)=4. 5 

340  THRESH (3)=2.34 

350  THRESH ( 4 ) : .8 1 

360  THRESH (5)=.09 

370  THRESH (6):. 09 

380  THRESH (7):. 81 

390  THRESH (8)=2.34 

400  THR  ESH (9  )  =  4. 5 

410  THRESH ( l 0 ) =7. 02 

420  THRESHU  1  )  =  7.02 

430  THRESH ( 12)=4.5 

440  THRESH (13)=2.34 

450  THRESH ( 1 4 ) = .8 l 

460  THRESH ( 1 5) = . 09 

470  THRESH ( 1 6 ) = . 09 

480  THR  ESH  ( 1  7 )  =  .8  1 

490  THRESH (18)=2. 34 

500  THRESH (19)=4.5 

510  THRESH (20) =7, 02 

520  GAMMA <1)=41 l .75 

530  GAMMA (2>=50S. 25 
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TRUCK  CONTINUED 


540  GAMMA (3) =58 7. 25 

550  GAMMA (4)=648. 

560  GAMMA <5)=668. 25 

570  GAMMA (6)=668.25 

580  GAMMA (7)=648. 

590  GAMMA <8)=587. 25 

600  GAMMA (9)=506. 25 

610  GAMMA (10)=41 1 .75 

620  GAMMA ( 1 1 )=41 1 .25 

630  GAMMA ( 1 2) =506. 25 

640  GAMMA (13)=587. 25 

650  GAMMA  < 1 4)  =  648 . 

660  GAMMA (15)= 668. 25 

670  GAMMA < 16)= 668. 25 

680  GAMMA  (17)  =  648. 

690  GAMMA ( 1 8 ) =587. 25 

700  GAMMA (19)=506. 25 

710  GAMMA (20)=41 1.25 

720  ZETA l=-5. 189 

730  ZETA2=0. 

740  ETA  1  =  . 028 

750  ETA  2=0. 

760  NU  l  =  - 1 . 08 

770  NU2=0. 

780  MU  1  =  - 1 . 245 

790  MU  2=0. 

800  NSP  =  0 

810  POINT  1  =  0. 

820  POINT2=0. 

830  SDZ  2SQ  =  0. 

840  T=0. 

850  NS TOP  =  0 

860  NPL=  1  2 

870  JJ=1 

880  DO  10  1=1,47 

890  10  Y (I )=0. 

900C - - - DATA  READ  IN - 

910  PRINT,"************TYPE  IN  THE  DATA  FOR  THE  FOLLOWING" 

920  PRINT, 

930  PRINT, "EXTERNAL  INPUT  DATA" 

940  PRINT, 

950  PRINT, "NAME  OF  PROFILE  INPUT  FILE" 

960  READ  2, FI  NAME 

970  2  FORMAT (A 6) 

980  PRINT, "DATA  VARIABLES’ 

990  PRINT, 

1000  PRINT, "TRUCK  VELOCITY  IN  INCHES  PER  SECOND  [REAL)" 

1010  READ ,VEL 

1020  PRINT, "PROGRAM  VARIABLES" 

1030  PRINT, 


TRUCK  CONTINUED 


1040  PRINT /NUMBER  OF  SPACES  BETWEEN  PROFILE  POINTS  [INTEGER] 
1050  READ  ,NSPACE 

1060  PRINT/NUMBER  OF  STEPS  IN  RKG  [INTEGER]" 

1070  READ  ,NSTEPS 

1080  PRINT/EXTERNAL  OUTPUT  DATA" 

1090  PRINT, 

1100  PRINT/NAME  OF  OUTPUT  DATA  FILE" 

1110  READ  2.FNAME2 

1120  PRINT, 

1130  PRINT ,"********£ND  OF  DATA  INPUT*********" 

1140  PRINT  3 

1150  XMPH  :VEL/1 7,6 

1160  3  FORMAT (///// , "T " , 1  OX , "Y ( 1 ) " , 7X , "AZETA  2" ,5X , 

11704  "ARMS",/) 

1180  WRITE(2;6) 

1190  6  FORMAT  (//,37X  ,  45  (  1H*)  ,/,37X  ,  1H*,43X  ,  1H*,/  ,37X  , 

12004  45H*  MURPHY  'S  M-37  TRUCK  PROGRAM  OUTPUT  FILE  *,/, 

1 2 1  0 &  3  7X,  1H*,43X  ,  1H*,/,37X  ,45(1H*)  ,///) 

1220  WRITE (2; t 1 )XMPH .NSTEPS 

1230  WRITE(2;9) 

1240  11  FORMAT (/,9H VELOCITY : ,F6.2, 15H  MILES  PER  HOUR.50X, 

12504  35HNUMBER  OF  STEPS  IN  RKG  I NTEGRATI ON: ,1 4) 

1260  9  FORMAT ( /  / , 13X ,9  ( 1H-) , 1 2HDI SP LAC EM ENT ,8( 1H-) ,3X , IOC 1H-) , 

12704  8HVEL0CI TY , 1 1 < 1H - ) , 3X ,9 ( 1H - ) , 1 2H ACCELER ATI  ON , 

12804  8C1H-),2X,7HC-G  RMS , / ,X , 4HTIME , 2X , 4HY ( 1 ) ,X , 

129  04  3(2X ,3HC -G , 4X ,5HPI TCH , 3X ,  5HFR-AX , 3X , 5HRE-AX  ,  2X  ) , 

13004  X , 5HACCEL , //) 

1310  DELTAT :3.07/VEL 
1320  CALL  OPENFd  .FINAME) 

1330  4  FORMAT (E20. 10) 

1340  50  POINT  1  :POI NT  2 

1350  IF(JJ-2)55, 40,55 

1360  55  READ ( 1 ,4  )POI NT  2 

1370  CALL  EOFTST (  1  , JJ  ) 

1380  GO  TO(60,40)JJ 
1390  40  POI NT 2=P0I NT  1 

1400  NSTOP :NST0P  +1 
1410  NSPACE:  1 

1420  60  PSTEP  =  (POI NT 2-POI NT  1 ) /NSPACE 

1430  NSP : 0 

1440  100  1:46 

1450  70  Y  (1+1  ):Y  (I  ) 

1460  1:1-1 

1470  I F (I )70,75,70 

1480  75  Y ( 1 ) :POI N T  1 

1490  P0INT1:P0INT1+PSTEP 

1500  NSP :NSP  +1 

1510  CALL  DIFFEQ 

1520  SDZ  2SQ :SDZ  2SQ+DZETA  2*DZETA  2 

1530  RMSDZ  2:SQRT ( (SDZ  2SQ*DELTAT )/T) 


TRUCK  CONTINUED 


1540  5  F0RMAT(6(X,G10.4)) 

1550  NPL=NPL+I 

1560  AZETA  2=DZETA  2/386. 

1570  ANU2=DNU2/386. 

1580  AMU2=DMU2/386. 

1590  RMSAZ 2=RMSDZ 2/386. 

1600  WRITE(2;7)T  .POINT  1  ,ZETA  1  ,ETA1  ,MU  1  ,NU  1  ,ZETA2,ETA2,MU2,NU2, 
161 0&  AZETA  2.DETA  2.AMU  2.ANU  2, RMSAZ  2 

1620  7  FORMAT (F6. 1  ,F5. 1 , 1 3F8. 3) 

1630  IF(NPL-54)120, 110,120 

1640  110  WRITE(2}9) 

1650  PRINT  5 ,T ,Y ( 1 ) .AZETA  2, RMSAZ  2 

1660  NPL=0 

1670  120  IF (NS TOP -47)80, 90, 80 
1680  80  T  =T  +DELTAT 

1690  IF(NSP-NSPACE)100,50,100 

1700  90  CALL  CLOSEF ( 2.FNAME2, 2) 

1702  PRINT  1 2, (I  BELL ,1  =  1,40) 

1704  12  FORMAT (40A  1) 

1710  PRINT  8.FNAME2 

1720  8  FORMATC///, "DETAILED  OUTPUT  IS  IN  FILE",X ,A6,/) 

1730  CALL  EXIT 

1740  END 
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DIFFEQ 


ICO 

SUBROUTINE  DIFFEQ 

110 

REAL  NU 1.NU2.MU 1 , MU  2, MASS .MASS  1  .MASS  2, 1 NRTIA 

120 

H=DELTAT/N STEPS 

130 

RH  =  1  ./H 

140 

INDEX=0 

150  100 

INDEX=INDEX+1 

160 

VAR  l=ZETA  1 

170 

VAR  2=ZETA  2 

180 

VAR  3=ETA l 

190 

VAR  4=ETA2 

200 

VAR  5  =  NU  1 

21  0 

VAR  6  =  NU  2 

220 

VAR  7=MU  1 

230 

VAR8=MU2 

240 

PZETA  2=ZETA  2 

25  0 

PETA  2=ETA  2 

260 

PNU2=NU2 

270 

PMU2=MU2 

280 

CALL  ALGEBR 

29  0 

FK  1  1=H*VAR2 

300 

FK 1 2=H*( 1 /MASS ) *(FORCK 1+C 1 1 *DSP DF 1 +FORCK  2+C  2 1 *DSPDF2-MASS*386. 0) 

310 

FK  1 3 =H *VAR  4 

320 

FK 1 4  =  H  *  C 1 /I NRTIA )*(A  *F0RCK 1+A*C 1 1*DSPDF1-B*F0RCK 

2-3*C  21 *DSPDF  2) 

330 

FK  15=H*VAR6 

340 

FK 16=H*( 1/MASS 1)*<-F0RCK l-C 1 1*DSPDF1+F0RCW1-MASS 

1*386.0) 

35  0 

FK  17  =  H*VAR8 

360 

FK 18=H* ( 1 /MASS  2)*(-F0RCK  2-C  21*DSPDF2+F0RCW2-MASS 

2*386.0) 

370 

VAR 1=ZETA 1+FK 1 1*.5 

380 

VAR  2=Z  ETA  2+FK  12*.  5 

39  0 

VAR  3  =  ETA  1+FK  13*.  5 

400 

VAR  4  =  ETA  2+FK  1  4*.  5 

410 

VAR5:NU1+FK 15*. 5 

420 

VAR  6  =  N U 2+FK  16*.  5 

430 

VAR  7:MU  1+FK  17*. 5 

440 

VAR  8=MU  2+FK 18*. 5 

450 

CALL  ALGEBR 

460 

FK  21=H*VAR  2 

470 

FK  22=H*(  1  /MASS  )*(F0RCK  1+C  1 1 *DSP DF 1 +F0RCK 2+C 2 1 *DSPDF2-MASS*38 6 . 0 ) 

48  0 

FK23=H*VAR4 

490 

FK  24-H*( 1 /I NRTIA )*<A*FORCK 1+A*C 1 1 *DSP DF 1 -B+FORCK 

2-B*C  21 *DSPDF  2) 

500 

FK25  =  H*VAR6 

510 

FK26=H*(  1/MASS  1 )*(-FORCK l-C 1 1 *DSPDF 1 +FORCW1 -MASS 

1*386.0) 

520 

FK27=H*VAR8 

530 

FK28=H*( 1/MASS2)*(-F0RCK2-C21*DSPDF2+F0RCW2-MASS 

2*386.0) 

540 

VAR  1=ZETA 1  +  . 29289322*FK  21  +  . 20710678 *FK 1 1 

550 

VAR  2=ZETA  2+. 29289322*FK  22+. 20710678 *FK 12 

560 

VAR3=ETA  1  +  . 29289322*FK  23+ . 2071 0678 *FK 13 

570 

VAR  4  =  ETA  2+. 29  289322*FK  24+. 207 1 0678*FK 1 4 

580 

VAR  5  =  NU 1  +  . 29289322*FK  25 +  .  207 i 0678*FK  1 5 

590 

VAR  6=NU2+. 29289322*FK  26+ . 207 i 06  78 *FK  16 

100 


DIFFEQ  CONTINUED 


600  VAR  7=MU 1+. 29289322*FK  27+. 207 1 Q678*FK 17 
610  VAR  8=MU2+. 29289322+FK  28+ . 207 1 0678*FK 18 

620  CALL  ALGEBR 

630  FK31=H*VAR2 

64  0  FK32=H*( 1/MASS )*(FORCK 1+C 1 1 *DSPDF 1 +FORCK  2+C  21 *DSPDF 2-MASS *38 6.0) 

650  FK33=H*VAR4 

660  FK  34=H*(  1/INRTIA  )*(A*F0RCK  1+A*C  1  1  *DSPDF  1  -B*FORCK  2-B*C 21  *DSPDF2> 

670  FK  35=H *VAR  6 

68  0  FK  36=H*( 1 /MASS  1 )*(-FORCK 1-C 1 1*DSPDFI+F0RCW1-MASS  1*38  6.0) 

690  FK37=H*VAR8 

700  FK38=H*( 1 /MASS  2)* (-F0RCK2-C 21 *DSPDF2+F0RCW2-MASS 2*386.0) 

710  VAR  UZETA 1-. 7071 0678*FK 21+1 .7071 0678*FK 31 
720  VAR  2=Z£TA  2- . 70  7 1 06  78*FK  22+ 1 . 70710678 *FK32 

730  VAR  3=ETA  1- . 70  71 06  78*FK 23+ 1 . 70 7 1  06  78+FK 33 

740  VAR4=ETA2-.70710678*FK24+1 . 707 1 0678*FK 34 

750  VAR5=NU 1-. 7071 0678*FK 25+1 .7071 0678*FK 35 

760  VAR6=NU 2-. 70710678 *FK  26+ 1.70710678*FK36 

770  VAR 7=MU  1-. 707 10678*FK 27+1 .7071 0678*FK 37 

78  0  VAR  8=MU  2- . 707 1  0678*FK28+1 . 70 7 1 0678*FK 38 

790  CALL  ALGEBR 

800  FK41=H*VAR2 

810  FK42=H*(  1/MASS  )*(FORCK  1+C  1  1  *DSPDF  1 +FORCK  2+C 2 1  *DSPDF2-MASS *38 6. 0 ) 

820  FK43=H*VAR4 

830  FK  4  4  =  H*<  1/INRTIA  )*(A*FORCK  1+A*C  1  1  *DSP  DF 1  -B  *FORCK  2-B*C  21 +DSPDF2) 

840  FK  45=H  *VAR  6 

85  0  FK  46=H*< 1 /MASS l)*(-FORCK 1 -C 1 1 *DSPDF 1 +FORC W 1 -MASS  1*386.0) 

860  FK47=H*VAR8 

8  70  FK48=H*< 1 /MASS  2) * ( -FORCK 2-C 21 *DSPDF2+F0RC W2-MASS 2*386. 0 ) 

880  ZETA lrZETA 1+FK  1  1  * .  1 6666  7+. 09 763 1 07*FK 21  +  . 569 03559*FK 3 1+FK 4 1 *. 1 668 
890  ZETA2:ZETA2+FK 1 2* . 1 66667+. 09 763 1 0 7*FK 22+. 5 6903559 *FK 32+FK 42*. 1 8 
900  ETA  UETA 1+FK 1 3*. 1 66667+. 09 763 1 07*FK 23+. 569 03559 *FK 33+FK 43* . 1 66S 

9  10  ETA  2  =  ETA  2+FK  14*. 16666 7+. 097 63 i07*FK  24+ . 56903559*FK  34+FK  44*.  1668 

920  NU 1 :NU 1+FK 1 5*. 1 66667+. 09 763 1 07*FK 25+ . 569 03559 *FK 35+FK 45*. 1 66667 

930  NU2=NU2+FK 16*. 1 66667+. 09 763 107*FK26+. 569 035 59*FK36+FK 46*. 166667 

940  MU1=MU 1+FK  17*. 16666 7+. 09 7631 07*FK  27+ . 569 03 559 *FK  37+FK  47*.  166667 

95  0  MU  2=MU  2+FK  18*. 1 66667+. 09 763 1 07*FK 28+ . 569  03559*FK 38+FK 48*. 1 66667 

960  DZ ETA  2= (ZETA  2-PZETA  2) *RH 

970  DETA2:(ETA  2-PETA  2) *RH 

980  DNU2=(NU  2-P  NU  2) *RH 

990  DMU 2= (MU 2-PMU  2)*RH 

1000  IF (INDEX-NSTEPS ) 1 00,200,200 

1010  200  RETURN 
1020  END 
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ALGEBR 


100  SUBROUTINE  ALGEBR 

110  REAL  NU1,NU2,MU1,MU2,MASS,MASS1,MASS2,INRTIA 

120  DO  100  1=1,10 

130  SEGDEF (I )=Y (I  )-VAR5-THRESH  (I > 

140  IF(SEGDEF(I ))50,100,100 

150  50  SEGDEF  (I  )  =  0. 

160  100  CONTINUE 

170  DO  200  I =1 1 ,20 

180  SEGDEF (I )=Y (I +26 ) -VAR 7-THRESH (I ) 

190  IF  (SEGDEF (I 1)150,200,200 

200  150  SEGDEF (I )=0. 

210  200  CONTINUE 
220  F0RCW1 =0. 

230  F0RCW2=0. 

240  DO  300  1=1,10 

250  F0RCW1=F0RCW1+GAMMA (I )*SEGDEF(I ) 

260  300  CONTINUE 
270  DO  400  1=11,20 

280  F0RCW2=F0RCW2+GAMMA (I  )*SEGDEF(I  ) 

290  400  CONTINUE 

300  SPDEF  1  =VAR  5 -VAR 1-A*SIN(VAR3) 

310  SPDEF 2=VAR  7-VAR 1+B+SI  N (VAR 3) 

320  DSPDF 1=VAR6-VAR2-A*VAR4*C0S(VAR3> 

330  DSPDF2=VAR8-VAR2+B*VAR4*C0S (VAR 3) 

340  FORCK 1=19.54*(SPDEF1**3)-192.42*SPDEF1*SPDEFI+913.55*SPDEF1 

35  0  FORCK  2= 1 . 39  * (SPDEF 2**3) -  1 . 24*SPDEF2*SPD£F2+307. 72*SPDEF2 

360  IFCDSPDFl >600,500,500 

370  500  C  1  lrCPOS  1 

380  GO  TO  700 

390  600  C 1 1=CNEG1 

400  700  IF(DSPDF2)800,900,900 

410  900  C 21 =CPOS 2 

420  RETURN 

430  800  C 2 1 =CNEG2 

440  RETURN 

450  END 
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